{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Content under Creative Commons Attribution license CC-BY 4.0, code under BSD 3-Clause License © 2018 by D. Koehn, notebook style sheet by L.A. Barba, N.C. Clementi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href=\"https://fonts.googleapis.com/css?family=Merriweather:300,300i,400,400i,700,700i,900,900i\" rel='stylesheet' >\n",
       "<link href=\"https://fonts.googleapis.com/css?family=Source+Sans+Pro:300,300i,400,400i,700,700i\" rel='stylesheet' >\n",
       "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' >\n",
       "<style>\n",
       "\n",
       "@font-face {\n",
       "    font-family: \"Computer Modern\";\n",
       "    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "}\n",
       "\n",
       "\n",
       "#notebook_panel { /* main background */\n",
       "    background: rgb(245,245,245);\n",
       "}\n",
       "\n",
       "div.cell { /* set cell width */\n",
       "    width: 800px;\n",
       "}\n",
       "\n",
       "div #notebook { /* centre the content */\n",
       "    background: #fff; /* white background for content */\n",
       "    width: 1000px;\n",
       "    margin: auto;\n",
       "    padding-left: 0em;\n",
       "}\n",
       "\n",
       "#notebook li { /* More space between bullet points */\n",
       "margin-top:0.5em;\n",
       "}\n",
       "\n",
       "/* draw border around running cells */\n",
       "div.cell.border-box-sizing.code_cell.running { \n",
       "    border: 1px solid #111;\n",
       "}\n",
       "\n",
       "/* Put a solid color box around each cell and its output, visually linking them*/\n",
       "div.cell.code_cell {\n",
       "    background-color: rgb(256,256,256); \n",
       "    border-radius: 0px; \n",
       "    padding: 0.5em;\n",
       "    margin-left:1em;\n",
       "    margin-top: 1em;\n",
       "}\n",
       "\n",
       "\n",
       "div.text_cell_render{\n",
       "    font-family: 'Source Sans Pro', sans-serif;\n",
       "    line-height: 140%;\n",
       "    font-size: 110%;\n",
       "    width:680px;\n",
       "    margin-left:auto;\n",
       "    margin-right:auto;\n",
       "}\n",
       "\n",
       "/* Formatting for header cells */\n",
       ".text_cell_render h1 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-style:regular;\n",
       "    font-weight: bold;    \n",
       "    font-size: 250%;\n",
       "    line-height: 100%;\n",
       "    color: #004065;\n",
       "    margin-bottom: 1em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\t\n",
       ".text_cell_render h2 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-weight: bold; \n",
       "    font-size: 180%;\n",
       "    line-height: 100%;\n",
       "    color: #0096d6;\n",
       "    margin-bottom: 0.5em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\t\n",
       "\n",
       ".text_cell_render h3 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "\tfont-size: 150%;\n",
       "    margin-top:12px;\n",
       "    margin-bottom: 3px;\n",
       "    font-style: regular;\n",
       "    color: #008367;\n",
       "}\n",
       "\n",
       ".text_cell_render h4 {    /*Use this for captions*/\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-weight: 300; \n",
       "    font-size: 100%;\n",
       "    line-height: 120%;\n",
       "    text-align: left;\n",
       "    width:500px;\n",
       "    margin-top: 1em;\n",
       "    margin-bottom: 2em;\n",
       "    margin-left: 80pt;\n",
       "    font-style: regular;\n",
       "}\n",
       "\n",
       ".text_cell_render h5 {  /*Use this for small titles*/\n",
       "    font-family: 'Source Sans Pro', sans-serif;\n",
       "    font-weight: regular;\n",
       "    font-size: 130%;\n",
       "    color: #e31937;\n",
       "    font-style: italic;\n",
       "    margin-bottom: .5em;\n",
       "    margin-top: 1em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h6 { /*use this for copyright note*/\n",
       "    font-family: 'Source Code Pro', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 9pt;\n",
       "    line-height: 100%;\n",
       "    color: grey;\n",
       "    margin-bottom: 1px;\n",
       "    margin-top: 1px;\n",
       "}\n",
       "\n",
       "    .CodeMirror{\n",
       "            font-family: \"Source Code Pro\";\n",
       "\t\t\tfont-size: 90%;\n",
       "    }\n",
       "/*    .prompt{\n",
       "        display: None;\n",
       "    }*/\n",
       "\t\n",
       "    \n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }  \n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"], \n",
       "                           equationNumbers: { autoNumber: \"AMS\", useLabelIds: true}\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Execute this cell to load the notebook's style sheet, then ignore it\n",
    "from IPython.core.display import HTML\n",
    "css_file = '../style/custom.css'\n",
    "HTML(open(css_file, \"r\").read())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2D SH finite difference modelling - the quarter plane problem\n",
    "\n",
    "Before modelling the propagation of Love-waves, let's take a look how accurate reflected body waves from the free-surface boundary can be modelled using the quarter plane problem. This is an important aspect, because one way to understand Love waves is by the interference of SH-body waves in a layered medium."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analytical solution of 2D SH quarter plane problem\n",
    "\n",
    "We seek an analytical solution to the 2D isotropic elastic SH problem: \n",
    "\n",
    "\\begin{align}\n",
    "\\rho \\frac{\\partial v_y}{\\partial t} &= \\frac{\\partial \\sigma_{yx}}{\\partial x} + \\frac{\\partial \\sigma_{yz}}{\\partial z} + f_y, \\notag\\\\\n",
    "\\frac{\\partial \\sigma_{yx}}{\\partial t} &= \\mu \\frac{\\partial v_{y}}{\\partial x},\\notag\\\\\n",
    "\\frac{\\partial \\sigma_{yz}}{\\partial t} &= \\mu \\frac{\\partial v_{y}}{\\partial z}. \\notag\\\\\n",
    "\\end{align}\n",
    "\n",
    "for wave propagation in a homogeneous medium, by setting density $\\rho$ and shear modulus $\\mu$ to constant values $\\rho_0,\\; \\mu_0$\n",
    "\n",
    "\\begin{align}\n",
    "\\rho(i,j) &= \\rho_0 \\notag \\\\\n",
    "\\mu(i,j) &= \\mu_0 = \\rho_0 V_{s0}^2\\notag\n",
    "\\end{align}\n",
    "\n",
    "at each spatial grid point $i = 0, 1, 2, ..., nx$; $j = 0, 1, 2, ..., nz$, in order to compare the numerical with the analytical solution. For a complete description of the problem we also have to define initial and boundary conditions. The **initial condition** is \n",
    "\n",
    "\\begin{equation}\n",
    "v_y(i,j,0) = \\sigma_{yx}(i+1/2,j,0) = \\sigma_{yz}(i,j+1/2,0) = 0, \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "so the modelling starts with zero particel velocity and shear stress amplitudes at each spatial grid point. As **boundary conditions**, we assume \n",
    "\n",
    "\\begin{align}\n",
    "v_y(0,j,n) &= \\sigma_{yx}(1/2,j,n) = \\sigma_{yz}(0,j+1/2,n) = 0, \\nonumber\\\\\n",
    "v_y(nx,j,n) &= \\sigma_{yx}(nx+1/2,j,n) = \\sigma_{yz}(nx,j+1/2,n) = 0, \\nonumber\\\\\n",
    "v_y(i,0,n) &= \\sigma_{yx}(i+1/2,0,n) = \\sigma_{yz}(i,1/2,n) = 0, \\nonumber\\\\\n",
    "v_y(i,nz,n) &= \\sigma_{yx}(i+1/2,nz,n) = \\sigma_{yz}(i,nz+1/2,n) = 0, \\nonumber\\\\\n",
    "\\end{align}\n",
    "\n",
    "for all time steps n. This **Dirichlet boundary condition**, leads to artifical boundary reflections. In the [previous notebook](http://nbviewer.jupyter.org/github/daniel-koehn/Theory-of-seismic-waves-II/blob/master/06_2D_SH_Love_wave_modelling/2_From_2D_acoustic_to_SH_FD_modelling_final.ipynb), we neglected these reflections. In this case, we want to incorporate them into the analytical solution, to test how accurate the boundary reflections actually are. We can describe the reflections by using **image points**\n",
    "\n",
    "<img src=\"images/SH_quarter_problem.jpg\" width=\"60%\">\n",
    "\n",
    "where we excitate the Green's function solution for the homogenous medium:\n",
    "\n",
    "\\begin{equation}\n",
    "G_{2D}(x,z,t,xsrc,zsrc) = \\dfrac{1}{2\\pi \\rho_0 V_{s0}^2}\\dfrac{H\\biggl((t-t_s)-\\dfrac{|r|}{V_{s0}}\\biggr)}{\\sqrt{(t-t_s)^2-\\dfrac{r^2}{V_{s0}^2}}}, \\nonumber \n",
    "\\end{equation}\n",
    "\n",
    "where $H$ denotes the Heaviside function, $r = \\sqrt{(x-x_s)^2+(z-z_s)^2}$ the source-receiver distance (offset), at image point ... \n",
    "\n",
    "- 1 (-xsrc, zsrc) to describe the reflection from the left free surface boundary\n",
    "- 2 (-xsrc, -zsrc) to describe the reflection from the lower left corner\n",
    "- 3 (xsrc, -zsrc) to describe the reflection from the bottom free surface boundary\n",
    "\n",
    "Using the **superposition principle** of linear partial differential equations, we can describe the Green's function solution for the source and reflected wavefield by\n",
    "\n",
    "\\begin{align}\n",
    "G_{2D}^{total}(x,z,t,xsrc,zsrc) &= G_{2D}^{source}(x,z,t,xsrc,zsrc) + G_{2D}^{image1}(x,z,t,-xsrc,zsrc)\\nonumber\\\\  \n",
    "&+ G_{2D}^{image2}(x,z,t,-xsrc,-zsrc) + G_{2D}^{image3}(x,z,t,xsrc,-zsrc)\\nonumber\\\\\n",
    "\\end{align}\n",
    "\n",
    "For a given source wavelet S, we get the displacement wavefield:\n",
    "\n",
    "\\begin{equation}\n",
    "u_{y,analy}^{total}(x,z,t) = G_{2D}^{total} * S \\nonumber \n",
    "\\end{equation}\n",
    "\n",
    "Keep in mind that the stress-velocity code computes the **particle velocities** $\\mathbf{v_{y,analy}^{total}}$, while the analytical solution is expressed in terms of the **displacement** $\\mathbf{u_{y,analy}^{total}}$. Therefore, we have to take the first derivative of the analytical solution, before comparing the numerical with the analytical solution:\n",
    "\n",
    "\\begin{equation}\n",
    "v_{y,analy}^{total}(x,z,t) = \\frac{\\partial u_{y,analy}^{total}}{\\partial t} \\nonumber \n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Import Libraries \n",
    "# ----------------\n",
    "import numpy as np\n",
    "from numba import jit\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from pylab import rcParams\n",
    "\n",
    "# Ignore Warning Messages\n",
    "# -----------------------\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "code_folding": [],
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Definition of modelling parameters\n",
    "# ----------------------------------\n",
    "xmax = 500.0 # maximum spatial extension of the 1D model in x-direction (m)\n",
    "zmax = xmax  # maximum spatial extension of the 1D model in z-direction(m)\n",
    "\n",
    "tmax = 1.12   # maximum recording time of the seismogram (s)\n",
    "\n",
    "vs0  = 580.   # S-wave speed in medium (m/s)\n",
    "rho0 = 1000.  # Density in medium (kg/m^3)\n",
    "\n",
    "# acquisition geometry\n",
    "xr = 125.0 # x-receiver position (m)\n",
    "zr = xr    # z-receiver position (m)\n",
    "\n",
    "xsrc = 250.0 # x-source position (m)\n",
    "zsrc = xsrc  # z-source position (m)\n",
    "\n",
    "f0   = 40. # dominant frequency of the source (Hz)\n",
    "t0   = 4. / f0 # source time shift (s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To build the FD code, we first define the velocity update $v_y$ ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Particle velocity vy update\n",
    "# ---------------------------\n",
    "@jit(nopython=True) # use JIT for C-performance\n",
    "def update_vel(vy, syx, syz, dx, dz, dt, nx, nz, rho, op):\n",
    "    \n",
    "    # 2nd order FD operator\n",
    "    if (op==2):   \n",
    "        for i in range(1, nx - 1):\n",
    "            for j in range(1, nz - 1):\n",
    "            \n",
    "                # Calculate spatial derivatives  (2nd order operator)           \n",
    "                syx_x = (syx[i,j] - syx[i - 1,j]) / dx\n",
    "                syz_z = (syz[i,j] - syz[i,j - 1]) / dz\n",
    "            \n",
    "                # Update particle velocities\n",
    "                vy[i,j] = vy[i,j] + (dt/rho[i,j]) * (syx_x + syz_z)                               \n",
    "                \n",
    "    return vy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... update the shear stress components $\\sigma_{yx}$ and $\\sigma_{yz}$ ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Shear stress syx, syz updates\n",
    "# -----------------------------\n",
    "@jit(nopython=True) # use JIT for C-performance\n",
    "def update_stress(vy, syx, syz, dx, dz, dt, nx, nz, mux, muz, op):\n",
    "    \n",
    "    # 2nd order FD operator\n",
    "    if(op==2):\n",
    "        for i in range(1, nx - 1):\n",
    "            for j in range(1, nz - 1):\n",
    "            \n",
    "                # Calculate spatial derivatives (2nd order operator)\n",
    "                vy_x = (vy[i + 1,j] - vy[i,j]) / dx\n",
    "                vy_z = (vy[i,j + 1] - vy[i,j]) / dz\n",
    "            \n",
    "                # Update shear stresses\n",
    "                syx[i,j] = syx[i,j] + dt * mux[i,j] * vy_x\n",
    "                syz[i,j] = syz[i,j] + dt * muz[i,j] * vy_z                    \n",
    "    \n",
    "    return syx, syz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... and harmonically averaging the shear modulus ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Harmonic averages of shear modulus\n",
    "# ----------------------------------\n",
    "@jit(nopython=True) # use JIT for C-performance\n",
    "def shear_avg(mu, nx, nz, mux, muz):\n",
    "    \n",
    "    for i in range(1, nx - 1):\n",
    "        for j in range(1, nz - 1):\n",
    "            \n",
    "            # Calculate harmonic averages of shear moduli        \n",
    "            mux[i,j] = 2 / (1 / mu[i + 1,j] + 1 / mu[i,j])\n",
    "            muz[i,j] = 2 / (1 / mu[i,j + 1] + 1 / mu[i,j])\n",
    "            \n",
    "    return mux, muz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we define a function to calculate the Green's function at arbritary source positions or image points ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Green_2D(ir, jr, xsrc, zsrc, x, z, vs0, rho0, time, G, nt):\n",
    "    \n",
    "    # calculate source-receiver distance\n",
    "    #r = np.sqrt((x[ir] - x[isrc])**2 + (z[jr] - z[jsrc])**2)\n",
    "    r = np.sqrt((x[ir] - xsrc)**2 + (z[jr] - zsrc)**2)\n",
    "    \n",
    "    for it in range(nt): # Calculate Green's function (Heaviside function)\n",
    "        if (time[it] - r / vs0) >= 0:\n",
    "            G[it] = 1. / (2 * np.pi * rho0 * vs0**2) * (1. / np.sqrt(time[it]**2 - (r/vs0)**2))\n",
    "    \n",
    "    return G"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... and assemble the 2D SH FD code ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "code_folding": [
     42
    ]
   },
   "outputs": [],
   "source": [
    "# 2D SH Wave Propagation (Finite Difference Solution) \n",
    "# ---------------------------------------------------\n",
    "def FD_2D_SH_JIT(dt,dx,dz,f0,xsrc,zsrc,op):\n",
    "        \n",
    "    nx = (int)(xmax/dx) # number of grid points in x-direction\n",
    "    print('nx = ',nx)\n",
    "    \n",
    "    nz = (int)(zmax/dz) # number of grid points in x-direction\n",
    "    print('nz = ',nz)\n",
    "            \n",
    "    nt = (int)(tmax/dt) # maximum number of time steps            \n",
    "    print('nt = ',nt)\n",
    "    \n",
    "    ir = (int)(xr/dx)      # receiver location in grid in x-direction    \n",
    "    jr = (int)(zr/dz)      # receiver location in grid in z-direction\n",
    "    \n",
    "    # half FD-operator length\n",
    "    noph = (int)(op/2)\n",
    "\n",
    "    # Source time function (Gaussian)\n",
    "    # -------------------------------\n",
    "    src  = np.zeros(nt + 1)\n",
    "    time = np.linspace(0 * dt, nt * dt, nt)\n",
    "\n",
    "    # 1st derivative of a Gaussian\n",
    "    src  = -2. * (time - t0) * (f0 ** 2) * (np.exp(- (f0 ** 2) * (time - t0) ** 2))\n",
    "\n",
    "    # Analytical solution\n",
    "    # -------------------\n",
    "    G        = time * 0.\n",
    "    G1       = time * 0.\n",
    "    G2       = time * 0.\n",
    "    G3       = time * 0.\n",
    "    vy_analy = time * 0.\n",
    "\n",
    "    # Initialize coordinates\n",
    "    # ----------------------\n",
    "    x    = np.arange(nx)\n",
    "    x    = x * dx       # coordinates in x-direction (m)\n",
    "\n",
    "    z    = np.arange(nz)\n",
    "    z    = z * dz       # coordinates in z-direction (m)\n",
    "    \n",
    "    # calculate 2D Green's function for direct SH wave from source position    \n",
    "    isrc = (int)(xsrc/dx)  # source location in grid in x-direction\n",
    "    jsrc = (int)(zsrc/dz)  # source location in grid in x-direction\n",
    "    \n",
    "    # calculate source position from isrc and jsrc\n",
    "    xsrcd = isrc * dx\n",
    "    zsrcd = jsrc * dz\n",
    "    \n",
    "    G = Green_2D(ir, jr, xsrcd, zsrcd, x, z, vs0, rho0, time, G, nt)\n",
    "    \n",
    "    # calculate 2D Green's function for image points\n",
    "    # shift image points by half the FD operator size to compare with FD solution\n",
    "    G1 = Green_2D(ir, jr, -xsrcd + noph*dx, zsrcd, x, z, vs0, rho0, time, G1, nt)\n",
    "    G2 = Green_2D(ir, jr, xsrcd, -zsrcd+noph*dz, x, z, vs0, rho0, time, G2, nt)\n",
    "    G3 = Green_2D(ir, jr, -xsrcd+noph*dx, -zsrcd+noph*dz, x, z, vs0, rho0, time, G3, nt)\n",
    "    \n",
    "    G = G + G1 + G2 + G3 \n",
    "    \n",
    "    Gc   = np.convolve(G, src * dt)\n",
    "    Gc   = Gc[0:nt]\n",
    "    \n",
    "    # compute vy_analy from uy_analy\n",
    "    for i in range(1, nt - 1):\n",
    "        vy_analy[i] = (Gc[i+1] - Gc[i-1]) / (2.0 * dt)           \n",
    "    \n",
    "    # Initialize empty pressure arrays\n",
    "    # --------------------------------\n",
    "    vy    = np.zeros((nx,nz)) # particle velocity vy\n",
    "    syx   = np.zeros((nx,nz)) # shear stress syx\n",
    "    syz   = np.zeros((nx,nz)) # shear stress syz        \n",
    "\n",
    "    # Initialize model (assume homogeneous model)\n",
    "    # -------------------------------------------    \n",
    "    vs    = np.zeros((nx,nz))\n",
    "    vs    = vs + vs0      # initialize wave velocity in model\n",
    "    \n",
    "    rho   = np.zeros((nx,nz))\n",
    "    rho   = rho + rho0    # initialize wave velocity in model\n",
    "    \n",
    "    # calculate shear modulus\n",
    "    mu    = np.zeros((nx,nz))\n",
    "    mu    = rho * vs ** 2 \n",
    "    \n",
    "    # harmonic average of shear moduli\n",
    "    # --------------------------------\n",
    "    mux   = mu # initialize harmonic average mux \n",
    "    muz   = mu # initialize harmonic average muz\n",
    "\n",
    "    mux, muz = shear_avg(mu, nx, nz, mux, muz)\n",
    "    \n",
    "    # Initialize empty seismogram\n",
    "    # ---------------------------\n",
    "    seis = np.zeros(nt)     \n",
    "    \n",
    "    # Time looping\n",
    "    # ------------\n",
    "    for it in range(nt):\n",
    "    \n",
    "        # Update particle velocity vy\n",
    "        # ---------------------------\n",
    "        vy = update_vel(vy, syx, syz, dx, dz, dt, nx, nz, rho, op)\n",
    "\n",
    "        # Add Source Term at (isrc,jsrc)\n",
    "        # ------------------------------\n",
    "        # Absolute particle velocity w.r.t analytical solution\n",
    "        vy[isrc,jsrc] = vy[isrc,jsrc] + (dt * src[it] / (rho[isrc,jsrc] * dx * dz))\n",
    "        \n",
    "        # Update shear stress syx, syz\n",
    "        # ----------------------------\n",
    "        syx, syz = update_stress(vy, syx, syz, dx, dz, dt, nx, nz, mux, muz, op)                \n",
    "                           \n",
    "        # Output of Seismogram\n",
    "        # -----------------\n",
    "        seis[it] = vy[ir,jr]   \n",
    "        \n",
    "    # Compare FD Seismogram with analytical solution\n",
    "    # ---------------------------------------------- \n",
    "    # Define figure size\n",
    "    rcParams['figure.figsize'] = 12, 5\n",
    "    plt.plot(time, seis, 'b-',lw=3,label=\"FD solution\") # plot FD seismogram\n",
    "    Analy_seis = plt.plot(time,vy_analy,'r--',lw=3,label=\"Analytical solution\") # plot analytical solution\n",
    "    plt.xlim(time[0], time[-1])\n",
    "    plt.title('Seismogram')\n",
    "    plt.xlabel('Time (s)')\n",
    "    plt.ylabel('Amplitude')\n",
    "    plt.legend()\n",
    "    plt.grid()\n",
    "    plt.show()     \n",
    "    \n",
    "    return vy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we compare the analytical solution with the FD modelling result using the spatial 2nd order operator and a wavefield sampling of $N_\\lambda = 12$ gridpoints per minimum wavelength:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "nx =  827\n",
      "nz =  827\n",
      "nt =  1520\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtMAAAFNCAYAAADCcOOfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xd4VFX6wPHvmUmZhDRIQk2A0EILEJAOitiwrUhRlFVBV9x1/amLa99VFl11xV53rVgRFcECCigEpYi00JskAUJNQhIS0mfO7487mUlImyQzmZC8n+eZh3PLOfeQm/LOmXPfo7TWCCGEEEIIIWrP5O0OCCGEEEIIca6SYFoIIYQQQog6kmBaCCGEEEKIOpJgWgghhBBCiDqSYFoIIYQQQog6kmBaCCGEEEKIOpJgWgghGgmlVK5Sqou3+yGEEMJ1EkwLIYSbKaVGKaXWKqWylVKnlFJrlFKDa6qntQ7SWic1RB+FEEK4h4+3OyCEEE2JUioE+A74C/A54AeMBgq92S93UEqZtdZWb/dDCCEaExmZFkII9+oBoLWep7W2aq3ztdbLtNbbAJRStyqldiulMpVSS5VSnUorKqW0UqqbvXyFUmqXUipHKXVEKfV3+/4xSqlUpdQDSqmTSqljSqnx9vP32UfCHynTpr9S6iWl1FH76yWllH+Z4w/Y2ziqlPrTWX2Yq5R6Uym1RCl1BrhQKXWlUmqLUuq0UuqwUmpWmbY62+tPtx/LVEr9WSk1WCm1TSmVpZR6zbNffiGEaFgSTAshhHvtA6xKqQ+UUpcrpVqWHlBKjQceASYAkcAvwLwq2nkXuENrHQz0BVaUOdYWsAAdgMeAt4E/AoMwRsEfKzP3+lFgGDAA6A8MAf5h7884YCZwMdANuKCSftwI/BsIBlYDZ4CbgTDgSuAv9v9XWUOB7sD1wEv2PlwM9AGuU0pVdh0hhDgnnXPBtFLqPftozA43tHWhUiqxzKugkj8KQgjhMq31aWAUoDGC3DSl1DdKqTbAHcDTWuvdWusS4ClgQNnR6TKKgd5KqRCtdabWevNZx/6ttS4GPgMigJe11jla653ATqCf/dypwGyt9UmtdRrwL+Am+7HrgPe11ju11nn2Y2f7Wmu9Rmtt01oXaK0TtNbb7dvbMN4MnB0cP2E/dxlG8D3Pfv0jGG8g4l37agohRON3zgXTwFxgnDsa0lqv1FoP0FoPAMYCecAyd7QthGi+7MHyNK11FMaocnuMEdpOwMv26Q5ZwClAYYwwn20icAVwUCm1Sik1vMyxjDJzl/Pt/54oczwfCLKX2wMHyxw7aN9XeuxwmWNly5XuU0oNVUqtVEqlKaWygT9jBPNlnd2XqvomhBDnvHMumNZa/4zxB8hBKdVVKfWDUmqTUuoXpVTPOjQ9CfjePjojhBBuobXegzEI0BcjML1Dax1W5hWgtV5bSb0NWutrgNbAIoyHGeviKEYQX6qjfR/AMSCqzLHoyv4LZ21/CnwDRGutQ4H/YrwhEEKIZumcC6ar8Bbwf1rrQcDfgTfq0MYUqp67KIQQLlFK9VRK3aeUirJvRwM3AL9iBJ4PK6X62I+FKqUmV9KGn1JqqlIq1D6V4zRQ1ywa84B/KKUilVIRGHOsP7Yf+xyYrpTqpZQKtB+rSTBwSmtdoJQagjGnWgghmq1zPjWeUioIGAF8oZRjcMTffmwCMLuSake01peVaaMdEAcs9WxvhRDNQA7GA3gzlVJhQBZGqrz7tdan7b+zPrPPk84GlgNfVNLOTcBrSikzsBfjAcO6eBIIAbbZt7+w70Nr/b1S6hVgJWADnrBft7o0fncCz9uzcqzCCMjD6tg3IYQ45ymtz/4Er/FTSnUGvtNa97XndN2rtW5Xj/buAfporWe4qYtCCHHOUUr1AnYA/vYHJIUQQtTgnJ/mYX9yPrn0o1Jl6F/LZm5ApngIIZohpdS19mklLYH/AN9KIC2EEK4754JppdQ8YB0Qa1+44DaM1E+3KaW2YqSEuqYW7XXGeOhmlft7K4QQjd4dQBpwAGNe9l+82x0hhDi3nJPTPIQQQgghhGgMzrmRaSGEEEIIIRoLCaaFEEIIIYSoo3MqNV5YWJju1q2bt7shPOTMmTO0aNHC290QHiD3tmmT+9t0yb1t2uT+Vm/Tpk3pWuvIms47p4LpNm3asHHjRm93Q3hIQkICY8aM8XY3hAfIvW3a5P42XXJvmza5v9VTSh105TyZ5iGEEEIIIUQdSTAthBBCCCFEHUkwLYQQQgghRB2dU3OmhRBCCCHcrbi4mNTUVAoKCrzdlQYVGhrK7t27vd0Nr7NYLERFReHr61un+hJMCyGEEKJZS01NJTg4mM6dO6OU8nZ3GkxOTg7BwcHe7oZXaa3JyMggNTWVmJiYOrUh0zyEEEII0awVFBQQHh7erAJpYVBKER4eXq9PJSSYFkIIIUSzJ4F081Xfey/BtBBCCCGEl5nNZgYMGOB4paSkkJCQQGhoKPHx8cTGxnL++efz3Xff1ftaKSkp9O3bt8bznnrqqXLbI0aMqPe1myKZMy2EEEKcIwoKIDsb2rTxdk+EuwUEBJCYmFhuX0pKCqNHj3YE0ImJiYwfP56AgAAuuugij/fpqaee4pFHHnFsr1271uPXPBfJyLQQQghxDvj+e4iOhrZtYdo0sFq93SPR0AYMGMBjjz3Ga6+9VuHYqlWrHKPa8fHx5OTkoLXm/vvvp2/fvsTFxTF//vwK9ebOnctdd93l2L7qqqtISEjgoYceIj8/nwEDBjB16lQAgoKCAKpst3RFxUmTJtGzZ0+mTp2K1toTX4pGxWsj00opC/Az4G/vx5da68e91R8hhBCisUpPhylT4PRpY/uDD2D4cLjjDu/2S7hPaeAKEBMTw8KFCys9b+DAgcyZM6fC/ueee47XX3+dkSNHkpubi8Vi4auvviIxMZGtW7eSnp7O4MGDOf/8813qzzPPPMNrr71WYbQcqLbdLVu2sHPnTtq3b8/IkSNZs2YNo0aNcvXLcE7y5sh0ITBWa90fGACMU0oN82J/hBBCiEbpjTecgXSpl1+GZjDo1+CU8tyrOqXTPBITE6sMpIEqR3pHjhzJzJkzeeWVV8jKysLHx4fVq1dzww03YDabadOmDRdccAEbNmyoz5cHoNp2hwwZQlRUFCaTyTH3u6nzWjCtDbn2TV/7S34tCCGEEGf5+OOK+3bvNl6iedmyZQu9evWqsP+hhx7inXfeIT8/n2HDhrFnzx6Xplj4+Phgs9kc266kiKuuXX9/f0fZbDZTUlJSY3vnOq8+gKiUMgObgG7A61rr9ZWcMwOYARAZGUlCQkKD9lE0nNzcXLm/TZTc26ZN7q9nHT1qYf9+44Nbi8VK//5ZrF8fDsD//refa6894rFrN5d7GxoaSk5Ojn3Lc4uYOK/h2vG8vDxKSkoc+3fs2MHs2bN59dVXK5yblJREly5duPPOO/nll1/YsmULgwcP5r333mPChAlkZmayatUqHn/8cXJzc7HZbFitVlq3bs2mTZvIzs7m6NGj/Pbbb+Tl5ZGTk4Ovry+nTp0qtzJgTk5Ole3u27evXH+LioooKCio8f/dGBQUFNT5e92rwbTW2goMUEqFAQuVUn211jvOOuct4C2A2NhYPWbMmIbvqGgQpQ8uiKZH7m3TJvfXs/73P2d57FgzV1wRznr70NPJk90ZM6a7x67dXO7t7t27HSsBenbqTPWB+tmrEQYGBrJu3TrOP/988vLyaN26Na+++ipXX311hbrvvPMOK1euxGw207t3byZMmICfnx+JiYmMGjUKpRRz5syhW7dupKSkYDKZMJvNXHLJJXzwwQeMGDGCvn37MnDgQAIDAwkODmbGjBmMHDmSgQMH8sknnzj6eOONN1babmpqKj4+Po7/h5+fHxaL5ZxYZdFisRAfH1+nuqqxPGWplHocOKO1fq6qc2JjY/XevXsbsFeiITWXX9rNkdzbpk3ur2dNmgQLFhjll16CoUONhw8B+vSBHTuqrltfzeXe7t69u9KpE02dLCfuVNn3gFJqk9b6vJrqem3OtFIq0j4ijVIqALgY2OOt/gghhBCNjdawZo1z+6KLoOxaG3v2QGFhw/dLCOHkzWwe7YCVSqltwAZguda6/sv6CCGEEE3EkSNw/LhRDgqCXr0gyL+Ynl2KACPXtDyEKIR3eTObxzatdbzWup/Wuq/Wera3+iKEEEI0RmWzmA0aBOY5z0B4OBsOt+FW3gVg2zYvdU4IAcgKiEIIIUSjVTaYvqnVYnj4YcjJIag4i7eYQT+2ysi0EF4mwbQQQgjRSG3Z4ixfvfs/5Y6ZsXEfz9MM1sQQolHzamo8IYQQQlStNFNHL3bRes8vjv0FbTsz8/j9fM51dEv2UueEEICMTAshhBCNUnY2pKYa5ammz5wHJk7k5Ppk3uROMoggWYLpJmPhwoUopdizp37JzaZNm8aXX35Z7TlPPfVUue0RI0bU6VqzZs3iueeqzGrskoSEBK666qpqz8nKyuKNN95wbB89epRJkybV67ruIsG0EEII0Qjt3OksX+2/zLlx/fV06AA+9s+WT56EvLyG7ZvwjHnz5jFq1Cg+++yzmk+up7OD6bVr13r8mvVxdjDdvn37Gt8wNBQJpoUQQohGqDSYDiOTvgX2JxGVgrFjMZuhY0fnuTJv+tyXm5vLmjVrePfdd8sF06UL50yaNImePXsydepUShfcmz17NoMHD6Zv377MmDGDsxfi++mnn7j22msd28uXL2fChAk89NBD5OfnM3LkSKZOnQpAUFCQ47xnn32WuLg4+vfvz0MPPQTA22+/zeDBg+nfvz8TJ04kr4Z3cF988QV9+/alf//+nH/++YCxZPf06dOJi4sjPj6elStXVqh39kh33759SUlJ4aGHHuLAgQMMGDCA+++/n5SUFPrak65X1e7cuXOZMGEC48aNo3v37jzwwAM13IW6kWBaCCGEaIRKg+kxJGDSNmPjvPMgPByAmBhjVxiZEkw3AYsWLWLcuHH06NGDVq1asXnzZsexLVu28NJLL7Fr1y6SkpJYY1/J56677mLDhg3s2LGD/Px8vvuu/HIdY8eOZffu3aSlpQHw/vvvM336dJ555hkCAgJYs2aNY5nwUt9//z2LFi1i/fr1bN261RGATpgwgQ0bNrB161Z69erFu+++W+3/Z/bs2SxdupStW7fyzTffAPD6668DsH37dubNm8ctt9xCQUGBS1+fZ555hq5du5KYmMicOXPKHauu3cTERObPn8/27duZP38+hw8fdul6tSHBtBBCCNEIlQbTbTlOiX+gsVFmae9/J93AQTqSSStO7kxr+A42ZbNmGZ8CuPKaMaNi/Rkzyp8za1aNl5w3bx5TpkwBYMqUKcybN89xbMiQIURFRWEymRgwYAAp9ndPK1euZOjQocTFxbFixQp2lp0bBCiluOmmm/j444/Jyspi3bp1XH755dX248cff2T69OkEBhrfc61atQJgx44djB49mri4OD755JMK1zrbyJEjmTZtGm+//TZWqxWA1atXc9NNNwHQs2dPOnXqxL59+2r82tSkunYvuugiQkNDsVgs9O7dm4MHD9b7emeTbB5CCCFEI1SayeO//IWZW26ne9FOCAlxHG9fmEQ0xiibaesW4FIv9FK4Q0ZGBitWrGDHjh0opbBarSilePbZZwHw9/d3nGs2mykpKaGgoIA777yTjRs3Eh0dzaxZsyod5Z0+fTpXX301FouFyZMn4+NTfeintUYpVWH/tGnTWLRoEf3792fu3LkkJCRU285///tf1q9fz+LFixkwYACJiYkVpqFUxsfHB5vN5th2ZeS6unYr+9q5m4xMCyGEEI1MdrZzGXF/f+jSwwf693fO7QByOvdzlP2S6pf9QXjXl19+yc0338zBgwdJSUnh8OHDxMTEsHr16irrlAaZERER5ObmVvkwXvv27Wnfvj1PPvkk06ZNc+z39fWluLi4wvmXXnop7733nmNO9KlTpwDIycmhXbt2FBcXV5gaUpkDBw4wdOhQZs+eTUREBIcPH+b888931N23bx+HDh0iNja2XL3OnTs7prhs3ryZZHu6muDgYHJyciq9livtepIE00IIIUQjs3+/s9ytG5jNFc+xde3uKAceO9AAvWpGZs0CrV17vfVWxfpvvVX+nBqmecybN6/cg4IAEydO5NNPP62yTlhYGLfffjtxcXGMHz+ewYMHV3nu1KlTiY6Opnfv3o59M2bMYPjw4Y4HEEuNGzeOP/zhD5x33nkMGDDA8TDgE088wdChQ7nkkkvo2bNntf8fgPvvv5+4uDj69u3L+eefT//+/bnzzjuxWq3ExcVx/fXXM3fu3HIjx6X/71OnTjFgwADefPNNevToAUB4eDgjR46kb9++3H///eXquNKuJylXhtwbi9jYWL13715vd0N4SOkTy6LpkXvbtMn9db9PP4XSGGf8eFi4sOI5+55eQI9HjDy7q0OuYFT2Yrf3o7nc2927d9OrVy9vd8Nj7rrrLuLj47ntttvK7c/JySE4ONhLvWpcKvseUEpt0lqfV1NdmTMthBBCNDKlI9NjWMm1OhP29oGuXZ3JpYGg/l0d5bZnZGRaVG7QoEG0aNGC559/3ttdabJkmocQwnXHjlFw9WTyuvej5O8PgQce5BBCQGmCg3t4mZu/ngg9e8L8+eXOCR/iDKajrclYi6wN2UVxjti0aRM///xzg057aG5kZFoI4RqrlVMjrqRVyhZj+/ntZJ+B0Def8W6/hGiCSkeme7PLubNPn3Ln+EcEc1K1prU+iT9FnNx2hNbndUQI0bBkZFoI4ZLMVz92BtLArwxl9rpLvNgjIZomrY2RaX8K6Ip9+oZSUEl2gqMW5+h09maZ6lEf59IzZMK96nvvJZgWQtRMawpmOUeg1zKc4azjha0XsWGDF/slRBOUnm6kxotlL2bs+Xa7dIGAgArnZoY4R6Lz9h9pqC42ORaLhYyMDAmomyGtNRkZGVgsljq3IdM8hBA1ylu3lXbZRh7bHIKY1nYpHDeS+n/yCVSTkUkIUUuuTPEolRfWAU4Y5ZLDxzzcs6YrKiqK1NRUx7LbzUVBQUG9gsimwmKxEBUVVef6EkwLIWqU/PRnlP4pXxV6Da/ODWbcOGP7xx+91i0hmqTShw/7UGa55jL5gcv6bdjdPLD3Vo7QgceHhDKoAfrXFPn6+hJTZkGc5iIhIYH4+Hhvd+OcJ8G0EKJGQau+c5Tz/3AdF44yMnSVlMDOnZCepomIrLj8rBCi9mozMu3TtZPjrLR0z/ZLCFE5mTMthKhW8aFjdMoxRsgK8SP+voto0QKm9t7C88xkI4M4/jfJ6CGEu9RmZDoiwllOl2BaCK/wWjCtlIpWSq1USu1WSu1USt3jrb4IIap24O2fHOXNlhF07dcCgKuiEpnJiwxiM9bfNnmre0I0Ofv3gx+FdON3Y4dSRp7pSkRGOsvNbLqvEI2GN0emS4D7tNa9gGHAX5VSlb/1FkJ4zcH1J8gkDIATfS5C2WdztLogznFOqyM7vNE1IZocrY1gulwmj5gYCAys9PyICPChmGgO0SZpnSykJIQXeC2Y1lof01pvtpdzgN1AB2/1RwhRuee5j1acogsHsP1phmN/9KW9sGFE1u3z9kNBgbe6KESTcfQo5OVBARY+8LsdRoyA4cOrPD8yEg7SiUN04r/bRsAxyeghRENrFHOmlVKdgXhgvXd7IoQoS2vYtAlAkUwX+l3c2nGsS1wLkukCgBkbWb/u8U4nhWhCSh8+3E8P3hjwFqxZAx9/XOX5ERFwgjbOHSdOeLiHQoizeT2bh1IqCFgA3Ku1Pl3J8RnADIDIyEgSEhIatoOiweTm5sr9bWSOH7dw6tQwAFq0KOHw4dWkpjqPl7ToSdczxqprv779LRayKm1H7m3TJvfXfb79th1grHQYGnqchITq36RarYriMsF04tLlZOXmuq0/cm+bNrm/7uHVYFop5YsRSH+itf6qsnO01m8BbwHExsbqMWPGNFwHRYNKSEhA7m/jsmCBszxkiA8XXjim3PFvY5bBjsUAhGcWM7iK+yf3tmmT++s+ixc7y6NGtWXMmLY11vnM710oMspdgtsS4sZ7Ife2aZP76x7ezOahgHeB3VrrF7zVDyFE1Yo/+oxbmEtftjM4vuKDTaZuXZzl5KSG7JoQTVJpWjyA7t1dq5Mb6ByZLkiRaR5CNDRvzpkeCdwEjFVKJdpfV3ixP0KIswxaMYe5TGc7/RgXsKrC8YA+zmA68GRyQ3ZNiCZp/37ozj6+5SouXPogLFpUY538YOezDCVHJJgWoqF5bZqH1no1IEumCdFI6cIiOuVsd2x3Gl9xydmWg5zBdMRpGZkWoj6sVvj9dxhPIlexGD5YDOlXwvjx1dYrCG0Dh42yPnGyAXoqhCirUWTzEEI0PseW78CPYgBSVAwxg1pVOCdqWBTF9vfkkSXH0WfyGrSPQjQlBw9CcfFZy4hXsfJhWSXhzmkepnQZmRaioXk9m4cQonE69t0m2tvLByMG0rmSz5Ei2vpwj/9rHCmMIIkuLM/xI6JFg3ZTiCajNsuIl6UjndM8fE9JMC1EQ5NgWghRqZLfNjvKZ3oOqvQcpWBV7B1s22ZsJx2CiJqTDwghKlGaY7rcyHSfPjXWM7VzjkxbTss0DyEamkzzEEJUKuzAJkc5cFTlwTQYKx2XSpZnEIWos337wJcielAmpUfPnjXW8+sQiQ1FNiHk+4aAzebBXgohzibBtBCiAl1UTKfT2xzb0dcMrPLc6Ghn+cgRT/ZKiKZt3z7oxu/4Yk9D2bEjBAfXWC80whd/CgkjmwfG7weT/GkXoiHJT5wQooJjP+3CQiEAh1RHYgZHVHluVFRpSXMiJd/znROiidq376z50i5M8QBo2RJK8AUgM9MTPRNCVEeCaSFEBccWO+dLHwwfWO1AV7/89eylB2dowdT5f2iA3gnR9BQWGtk8apvJAyAszFmWYFqIhicPIAohKij+1Tlf+kxs1fOlASKiA+iB8eRUcI7M8xCiLg4cAK3rPjJdSoJpIRqeBNNCiAq+5hpWY2EQm/AfM7zac8P7dXCUIwpSPd01IZqk0rR4dRmZbtkSojlEF5LofvQU7OsLPXp4oJdCiMpIMC2EKEdreOfgJaRzCQD7bqn+/HZ9WpGPhQAKCNY52LJOYwoLaYCeCtF0lAbTt/Eu94zdwY0DdkGvXi7VbdkS7uN57uEVyAC+ex5mzvRcZ4UQ5cicaSFEOYcPQ3q6UQ4Jga5dqz8/IFBxzOQcnc7YKqPTQtRWaY7p3xhKxvjb4PnnjR9AFwQHQ6ZyrlBqTTvliS4KIaogwbQQopyNG53lgQNdy7KVZWnnKJ/addwDvRKiadtXJrV0bWdomEyQHxDu2C48luGmXgkhXCHBtBCinE3OZw8ZVP2zhw5ngpwrsJ1JlhXYhKitPXuc5bpMdy4McgbTJSckmG7utIajR6GgwNs9aR4kmBZClDP+rcv5mdG8yL2MinEtO0dhmDOYLk494amuCdEkpafDyZMQRA6BAZpOnWrfRkmoM5jWaRJMN2e7dkG/ftChA0REwKuvertHTZ8E00IIB11ipXf6L4xmNffyMnH9XfsVYY1wBtPWYzIyLURt7LRnw/uNIaQWt8Z08Vjj4YVa0K2cwbTKlGC6uTp0CC66CHbsMLbPnIG774Z587zbr6ZOgmkhhMOxVftowRmjrNoRM6JdDTUMqk1rR9mUJiPTQtTGrl3gTwHd2U/LknRISIDw8BrrlaUinOebsySYbq7++lc4XsljK/feK1M+PEmCaSGEw5FvnSsfprSsfuXDsnyjnCPTfpkSTAtRGzt3Qi9244PV2NG1KwQG1qoNc2tnMO2XI8F0c7R8OXz3nVFWCla+sp2bwpfQigxOnoTPPvNu/5oyCaaFEA5F65xPH+b0cPHpQ8A6egzxbKYDqTwSu8ATXROiydq5E+LY7twRF1frNgIigyjCFwDf4nzIz3dX98Q54rnnjH/9KGR95+sZc3c/Psy4koN0YhJf8MYb3u1fUybBtBDCIWS/M5i2jBjocr3wbi1JJJ6jdOBImp8nuiZEk7VrV/2D6dAwRQZlpoZkyOh0c7JjByxbBqD5mJsYnPy541gQZ/iEqZRs2EyqLAPgERJMCyEA4+HDzllbHNtRf3A9mG7tnDLNCZnlIYTLSjN51DeYDgmBnfQhkf7s6TAWrFY39lI0dqUZO6bwGZP5wnkgKAgAP4r5Dw+yeLEXOtcMSDAthADg6Mq9BOscAE6oNnQeFeVy3bLBdHo6lJS4u3dCNE2lmTzcEUxfwo/Ek8gTF/xEnfLriXNSURF88YXxEOtz/N154PbbYeNGbCYzYHx/HPhwjZd62bRJMC2EAODIog2OclL4YExm5XJdX18Ib2kjnHR66Z1kHM7zRBeFaHJ27YJWZNCBo8YOiwW6dat1O6GhznJ2tps6J84JP/4ImZkwjbnO76N27Ywl6WNjyZl0K0u5lJv5gHk7+2Gzebe/TZFXg2ml1HtKqZNKqR3e7IcQAorXOoPp3J6Da13/x4KRpBPJTvqS80uiO7smRJO1Y8dZo9J9+oDZXOt2QkKc5dOn3dAxcc6YP9/49zw2Onfefz8EBwMQMu9/3Nx6KR9xM6nZweza5YVONnHeHpmeC4zzch+EEEBg8k5necyQWtfPC4hwlHMPyMRpIVyRmAgDcaakpH//OrXTWILprVvh2Wfh7bdlhLwhFBTAokVG+XbeYdfc3+CPf4QZMxznKJNi5EhnnTUy08PtvBpMa61/Bk55sw9CCONZpfOLf6IXu7iZD+g8ZVit2ygIdeaaLjwsqyAKUROr1Qimw8lwpLVjkOspKcsKDYWe7GY67zHp0AvGZ/8N7MUXIT4eHnzQiOViY431Z4TnLF3qfPPUpQv0unkwfPQRtGhR7ryywfTatQ3YwWbCx9sdqIlSagYwAyAyMpIE+clssnJzc+X+eklyciC5eUPYQy/Swrsw/eQ69qfVro1MX4ujfHLHrnL3Uu5t0yb3t24OHQokL28I/+DfvN7qQb55ah5FkZEU1eFrmZXly8Vs5lXuhkw48to17Pep/594V++RGfhJAAAgAElEQVRtYmIoM2fGl9t34gRcfrmV559PpHfvnHr3RVT0yiu9AGMgY9iwg6xalVzpeX5+oUA8oEldkUhCgvGxgfzsukejD6a11m8BbwHExsbqMWPGeLdDwmMSEhKQ++sdyWV+/44c6c+FF46pdRvLu26BfUa5lVUzssy9lHvbtMn9rZtPP3WW44eFcN4dd9S5rcJCWECSY7t9ixZ0cMM9ceXeag0PPODcjosz0v2dOAEFBWbmzBnE1q3lp6KI+svLg/XrndsPPNCJ/v0rz+IyMF7jf/cMruI72qceozDmGP6d2srPrpt4e860EKIR2OB89pAhtZ8uDZRfzticLQtGCFGTzWWmSg90Pa17pfz94YxPmGPblpFVvwZrYcUK5+8QiwUWL4aff3ZmGElJgbvvbrDuNBuLF8OZM5pPuYEnWr9Kv/bpVZ4bEqoYaNlFe44BcOSr9VWeK2pPgmkhBEGL59OdfYCuczDt284ZTPvlSDAtRE22ONdIIj6+6vNcVdzCGUxbTzVcMF12hP3WWyE6Gnr0gP/9z7n/gw9g9eoG61KzMH8+DGITN/AZ/zh5Nyq2h5F0ugpHop3PwuT++GtDdLHZ8HZqvHnAOiBWKZWqlLrNm/0RojnK3n+SZw9NYR+xHKIjw4fUbeW0wChnMB2YV/UIiRACbDZjZPomPuRavmJQVP0z4NhCnMG0zmyYYLqwEBYscG7ffLOzfP31MHmyc/v//k8WZnSXnBxjZHoac507r74a/PyqrFM80BlMB2yTYNqdvJ3N4watdTutta/WOkpr/a43+yNEc7R/rjNP0qkWHQkKrX2OW4AWHZ3BdFChjEwLUZ09eyArS/M0D/MVE+k0tK2RdLoedKgzmFbZDRNM//qrMwVe5872aWLr18N118HQobynbmOYvzEEn5gIb73VIN1q8r79FnRBATdS5mOBadOqrRM41hlMRx3bIO9s3EimeQjRzOUv+8VRTu81us7thHRx5pkOLfFSMJ2ezunxN1MQ2JKCgDAyx02BAwe80xchqrF2LXTioHPFuuBg6NmzXm3qsJaOsjmnYYLplSud5UsvBfXRhzB8uLG+9W+/EfT5e6wpHszfmQPAP/9prNYn6mf+fLiab2mF/YvZuTNccEG1dTqOiOK4PfNHgPUMJCVVe75wnQTTQjRz4XucExkDLhlV53ZaxYRixUQRvuToIHRxiTu657q0NLJ6DiPk64+w5GdhKcim5dL5FPTsT8HC7xu2L0LUYO1aGEmZ1TOGDYN6prLzb9WCEoxPlswFedXOn3WXFSuc5WtjEuH22430HmWYbFbm8ACj+IWMDJg92+PdatKysuCHH86a4nHLLWCqPqTr1g12EOfYLtq0vZqzRW1IMC1EM3bm2Gl65DpTCnS/ZUSd27IEmmgfmI0/hXTgKKfzGjbz5uE/3ElYRsVRaEvJGUwTx5O3eGUltYTwjrVrYRRlnsgbVfc3sqVCQhVZOKd6eHoJwsLC8qnZLlz6kDOA79ULPv/csVrIzkmPsxrjk6/XXjOmuYi6+eorCC86yjh+cO685ZYa61kskBLaz7Gd9fM2T3SvWZJgWohmbPcbK/HBmDe30zKQyNhW9WrPEhEEKAAyGnCmR87StUT/+qVj+9EOc5k5dA0pGDlX/XQR+prxFB043HCdEqIK6emwd+9ZI9Nll6iro5AQygfTWZ6d6rFzpzN2vqb9b/gnLDU2TCbjqcTJk2H5cnj/fXrPf5zR9llkJSVw330e7VqT9umn8Ec+xozN2DFmDMTEuFQ3u6MzmC7aLCPT7iLBtBDNWP43yxzl43GX1ru9COe0adIbMKHH0fuec5QXBk7lgZ238MKvI1jx+M8coT2F+HG/9Wn+9nxUw3VKiCr88gu05gT9sAczPj4wdGi92w0NheVcwgImkDjwVmMo0oPKpvb7i//7zo0bbzRGpgECAmDaNJRJ8fLLoIz32ixZYkxVELVz7Bis+ElzG2XyNdTw4GFZOs4ZTAf+LiPT7iLBtBDNWMfdzmA6ZFL9g+lwZ0KPBhuZLimBO0/MYi63UIgfPPSQY7GIW2d1ZPGdSxjBWt7kTt54U/HRRw3TLyGqsmwZXIrzZ48RIyAoqN7thoTAnbzJJBbwydh3jYTPHlS66IyZEkaf+MJ5oIqVHOPjjTzUJqzcyes8eG8hxcUe7WKTM38+jGQ1saXLzQYHw6RJLtcPG96LEsxYMZFf4osqaeBnW5ooCaaFaKaOrkmmU/HvAJwhkL4z6j5fulSXFicYwnouZwmFu36vd3uuSEiAFen9mM5c+kce44oH+pY7fvtr/el23SDH9h13wNatDdI1ISq1bBlcxlLnjssuc0u7ZZfrPn3aLU1WqzSYtuLDhjc2wpNPGik9RlT9u+Spv6Wx3DyO17mLO/b+jddf93w/m5KPP4a2HHdk5eDGG6FFC5frd+trYQCJBJHLtT12oev50KswSDAtRDOV9PK3jvLO8AsICPOvd5uTU+awnmEs4UoiEhbUXMENvigzIHbplFb4n/XfUAreeceZdSw/H268Np/TiZIWSjS8AwcgOclWfmTaTcF06Scy4PlgWuvyabF7XNoZHn0Uli6tNqtE65+/ZKz1R8AYRd/14Af83jDvu895W7bApk3wJZPp6nuYnA++gnvuqVUbPXrATvpSQIB83d1Igmkhmql3Dl3KkzzKHmI5c/E1bmlTt3LO89ANMM9Da1i0yLlddrW1soKDjSfgg4KgK7/zcfIIckdehi3Ts9kOhDjb0qUQzxZak2bsiIhwz1rilB+Z9nAiD44ehdxcoxwWBm3buljxz3/GOuk6x+bLRX/mqclbsNnc38empuzy7OMn+xJ887XOuekuatfOOZU+MxNyc+u2SJcoT4JpIZqhtDT4eGNP/smT9GI3sXP+5JZ2TZHOYNqc6fknEHdvyCX9pJGNpFWr6hMi9OoFH7xVSAJjiCeR9nm/s3fUbRVy4grhSYsWwU76MJ6F7Bs4xZjvWkN+YFeFhMBQfuUJ/sFtW+4yJth6yN69znLPns4HC2ukFOb33yW/S28AAijgH4kTef7RU+7vZBOSkwOffOLcrmJaeo2UMtZ3KXX8eEC9+iUMEkwL0Qx9/bVzJdmRIxXto90zOuHb1hlM+572/Mj0mX8+TToRLOIa7uqzssaYZMIN/iRc5cz80WvXAnb/6blqagjhPmlpxiInhVj4mvEEfj0P3nzTbe0HB8MgNvEP/s3E468bDxR4SGme6EhOMqhTLd84BwURsOQrCvyCAehCMoOemcTSrwvc3Mum49NPwZqbBxgDA6PrvlgtXTrb6Ewyl/EDEd8vd1MPXXfyJMycCV26GN+zffvCI4/A8eMN3hW3kWBaiGaobEaLWjwIXiNLB2cwbTnj+WA6aOMqWpLFNXzDsFjX1iiesnAKi9rf6dju9d4D/P5vz43gCVFqwQLnm9hRoyDKzZkag4MbLs90aTB9N6/w2vxI6NqVWqXKiY3F99MPHZtjWQmTJrF9k+dXbTzXWK3wwvOaFYxlGZfw79E/oKj7J2q9o06TTBd+4HKu/HoWDTnHJjHRmNX04ouQnGxMFdq5E55+2vgWeuopI0PTuUaCaSGamf0rDnPsZyOtko8PXH+9+9oO7OgMplsUeDaY1sUldDrlXL0x5o+uLXrh4wPD173ABn/n+dH/uJlDH8gKicKzPvvMWXbnz12p4GDIpsxTiB6cOF06zWMIvxmFpCTw86tVG+aJ48n9xzOO7ctKFnNo5BT2b5cR6rIWLoSo/SsYxnou4UfGz73G+Jijjtr2DOMErQHwtRbBoUPu6mq1jh+HK6805tsDBHOa8SzkXl5kOu/RJW87jz6qGTWKc+7hSAmmhWhmjt33HPuIZQUX8sCI1bRr5762gzs6V1AMKnZtpLiuTiTsJpB8AA6raHqMbuNy3TYd/QlN+IZ9ZiPFhz9FRE67kqQ3l9ZQU4i62bULNq7KZTrvEajymTjR/dcIDobTNEx+vKQkAM1gNjh3DhlS63aCnniQ47c96ti+snAhy4Y9xqZN9e9jU2CzwZNPaJ7hIcc+NW0atG5d5zZjYmAfPZw79u2rRw9d95e/lAbSmocsL5Ee2JGFTOBFZvIet7GdfuyiNwPXv8GQAUV8+OG580iLBNNCNCPZB7Pon/gBABeSwPhx7h0BCunU0lEOtWV69Dfhka83OsrJ4efV+hmuHsNacebLHziq2gMQQD5Rd17Nzgc/rKGmELX36qtwMx/yHrdx3DeKdl+84vZrBARAjnKOTOtszwTTNpsxmNmRQ7TEPpWkVavyT7bVQtu3nyB1irG++Fb68XDePzj//PIj+c3Vhx9Cz23zGYzx+05bLPDPf9arzZgY2Eusc0fZp0k9ZMWK0sxLmv/yZ54u+Bt+eRU/OenFHt7gr2w404slt3zGH6fqBsmZXl8STAvRjCROf5lQjF9gSb49GHT/WLe2H9ImwFiFEGO0tzjHcx/XFv/qHLo603NQNWdWLX58J1I/TuCQ6giAH8X0efYWPr51xTk5b080ThkZ8OkHxdzLSwAEF52qRfoL1ykF1hbOkWlblmemeZw4AUVF0Jcyiabj4ur+f1KKqE/nsP9vb3Bj6BJyCCEvD264AW66yXhgrTnKyYH/PJTJS9zr2Kf+7//qPdn+7JFpvcfzwfS//23825M9TDeXGbDo3h3++le49lrjoxW7riTxJn9hybws4uNh/XqPd7FeJJgWopnITj7FgIQXHdtHpv8Tk497fwWYzIos5RydPn3Qc1M9gpKcyxgGjhxY53aG3NidnCWr2etjpOqaxxRuev9Chg6FDRtqqCyEC558Em7If5ce7AdAh4bCLbd45Fq2YOfItPLQkF5KivFvHNudO/v2rfRclylF9xf+woJfO9Cli3P3Jx/bWN9hAsuufZMzGc1rLvW998KDJ/5GW04AYGvX3lgYp57CwuBoYHfHduEezy5gtW2bMTINsN/ci4wFq6BNG2P1xh074LXXjIUAjh6F555DtzSmC/6Lx8miJUlJxgO7Tz1lvIlrjCSYFqKZ2DbuAUK1c1R6+MtTPHKdQ/7d2UMs6xhGdnqxR64B0D57t7N8SZ96tdVnXDShu35lXtT93M7bgGLzZmMK6LXXwo7vUs6dyXuiUdm/Hz5+LYtZzHLsUw8/XH6FFTdSIc7RPZV72iPftwcPGv+WG5mubzBt17OnsdLfTTcZ29OYy9UlC7l00Z3kRMbw46hZpKxo+quXfv45qPfeYRofOPaZXn+t/DKX9VAc7XzHYvvds1/Pd95xlidOhHbXDIGNG+H998s/tBoUBPfdh0o6AE89xYiP7nT8mJSUGO8jHug0n6XvH210v44lmBaiGdj09DJG73vXsX1i5rP4WHw8cq2/9PmFXuxhBOtIC+zkkWvkJqfRymZkC8mlBZ1HR9e7zbbdg5mc/CyP/Duo3JLkqxel0fnqviQH9mHt5BdI25BS72uJ5qG4GG6+SfNqyZ8do4s6Kgruvttj17QE+3KGQACUzQZnzrj9GqXBdLmR6bg4t7UfEmLMFf7uW80Dfi879rfVx7l4zb/ofFFXtgSPJuGKZ9nz+TZs1kYWWdXTmjUw74+LeQNnCk/++Efjnb2bmLrGOMr+R5M9lh6vpKT82kEzZtgLUVFVZ38JC4OHH+a6P/qRmAjDhhm749nMC8dv4MJbO/NN+HSWPPYr+XmN495LMC1EE5eybB9dHnWOQq9tP4nhz7hn+fDKtHTO8iDTQ7M8UpftcpRTLD3x9XPP/FMfH2PxgB07cGRbeIBnCeIMMQW7GfHlfUQOiWF/YH9+GfEAmx//muzf656iSjRdWsPf/w5D17/MFJzRhHr5ZeNJQQ9piIweKSmgsNHdPm0FgN693X6dK69SdD2xltUTX+SYuUO5Y/G5qxnz/YP0vL4/J/068GuHiXz0119ZvhxOncOLKX7/PYy/NI+3i2/BD+OTvZLeceXXEneDyK4hpBEBgLmkyJmvzs1WrHDOeW/XDsaMqV39mBj4+Wd4bo7mFfNMTGj8KOaazLlc8cRwjgd145u4R1jz9M+cTvfeHBDPDE0JIRqF/V9tp8V1V9BSG1HtcVM7uix+1aPXbIhgOm3zYXqgMKHJaN3L7e136wZffmnM9Uu9zkrO3iCCyXUc756/je7rtsG6OTAbDvt05kR4H/aNmAaTJtGtG3TqBJGRblspWpxDSkrggfs1plde4kVmOg/86U8wYYJHrx0cDO/wJywUcOWUEHp7IHA/eBDaccyRmpLw8PI/+G7kG9aCUV/eizXvL6x//BtMH31A/Ikf8MHqOKet7Rhtj37Fw2/cRcIb9n1toUcP+Mepv+HTKYqA2I4E9YqmVf9oIuPa4mtxz6qv7pKZCbNnw8svg9aB/Jn/8iWTKY6OwfeHxRAY6NbrRUdDEl2IxL56ZVKS+1cRApYsMVbJXMR4Dva6FXP+FGM6Ry34+sJ998HpyHtIeriYLsfWOo7F6CRidjwNjzxNziNB/BY6jNPdBpH7hxvpeFU/evcGi8Xd/6uKvBpMK6XGAS8DZuAdrfUzNVQRQrigqMhYpXjLgzuZa00FII8AMt77hj4D2nr02g0RTC9u+UcuZSI92MefrvDlAs9chn79oN+eF8hI+Rc/PTyfFsu+Iv7UT/hTfgQkuiSF6BMpfLpwLC8udO7394fPfafS0ecI+SFtKQ6LxBYeiblNBD7tIwmIjiQouiUBbUJo0S6EoKgwfPwb1x95UTubNhkPjo1dPZt/lZknrUeMQL30ksevHxwMj/EEAG2vgN4eiHEPHoRI0kilA1EcMZau8zBzoD9D50yGOZPJ2nuC7c9+j/nHH+h1eBktdSY2FBs5z3H+8eNQfDydS3gJdgCLnW2VYOaIqT3ZltbkBUZQGBRBcVgEP1/zAq0iTISHG5n+Wvrm0up0CpbIYALbBNOibTD+Qb5u+z9lZ8OmH9LY987PvPbreezMdU6L29BxEsenvkTbGdcYka+bRUfDZowHt4s6dGF0LQNcV/3wA/yRjxnBOkasWAcTPodly2rfkFKE3HItIbdcy+ml60h69B26bFlAiM2ZsSaYXIZk/wibfuS6TQP54vF+KGWMiHfuDI+n34V/eDC0a4c5qi2WTm0JiokkrHMYrWJC8QsNqHNGGq8F00opM/A6cAmQCmxQSn2jtd5VfU0hRGVO7c9g94ZcFm7uxLx5pZ/aXc9EPuECVrH/PwsZdMt5NTVTb931PqawiZZkErQlDhjt9mvs2gUFBLCN/kR6KpIuI7xzMBfN+xPwJzIP5bDh1Z8o+mk14fvW0uvMRsfHsTsp/yBkYSH0LVxHF5IhC6hhobHR/MxGy2hCQoygKDgY3kq5FJMZSnwDsPlasPpZsPkFoP0txssSAAEWTAEWDlxwK6bQYHx9jdEcP2s+7Xcsw+Tvi/L1wWzxxeTvi9nfx/jX4ovJz9jvE+ALnTvj62tMdzGbwWQrwVSQh8nHhMnHhNnP7Cgrswllcn96t3OJzQYnT2iSt2SxNjGQhUv8WbPGOJbHVY5g2jZkKKYlS6BFC4/3qWxMlJPj/va1NoLpMwwgmlQyUvNpRcPOqwiLbcPod6cB07AVlbDvm10cWrqbqeZgNm0yPlEqKjLSsFXGBysdbIfpkHcY8oB0OE0wFyaWf7MzlvX8xMXl9hXgT64KJs8UTLHZQonZn2KzhcNBvXi5/3v4+xtTgf38oHf2OkYmf2SsB261oa02KC7GN+cUAfkZRBUlMZaTjAWSeYadPAjAZZcZz+a1bXeP+794dh07wg28CcCgtrCx7gmRqpScbKSwfp8vnDvdsPRnyGXDGXDZcHT+6xx4/XvSP/qeqD3L6VCU4jhnE0a6VK2Nv4VpR4u4lNerbbcYH06bwsjzDeWeYb/h26ZVteeX5c2R6SHA71rrJACl1GfANUCVwbTt4ClWDfxbpcfOfrLzu5FPU+LjHNsPOnOCSzY8VWVnylYvMfvzzchnyx2PyNzPBVtdS7KfE9CaHwaXT6oefXwDw/bOrfniwMnQ7qzs78wrqTX0TP2R+KQFrlTnYORg1va8tVz9gUlf0it1uUv9393hYjZ1mVyu/ug9b9M57beaLw5sjJnMrqhLHXUBxm17lrbZ1eeyLMzL45fAj1gVO4OkyKHl6k/a9DAh+cdduv6SuAc5FtqzXP1p6+7A1+paWqX5A/9DdkBbR13fknxuXT+j4olVPPfw7tC3KPYJcNQPyzvKjVvud+naRWYLbw99t1zfO2TvYsLOJyqcay7Kx5J3isDCU4QVnSRSp3GS8TxPmaFRFC/1eJPYf59i0KR+LvWhvgYe/Y6/YyzA8PPme/BUMF2ql/tneVSrZcdgRs0ZD4wHIC+zkAPL95H+8y5G+Q4mINmYU3roEGRnWoki1eW2cwimoAAKCox5hiasDMa1n1uA6xdezwmc2RyiSeOQvZ81sWIq99E5wDA2sI4RVdaxobBhYgQm9qvODAjYh8mE4zWu5DteOXMrWpmwYkYrEzZMaExYlbnMKJBil2Ug93X4DKWcu6/K/oQ/pT8NGDu0UmgUoMqUjfq/BV/Eq9HPOuoqBRNPvskV6R+ilSrXRtn6CtAofo6YwFdRdzvqKgXXH5rDeZnL0UphthZjshZhthbhYy3CbCvCv+QMbTjBcAqZzRLWcLnja5NoGsSOPjfQe1RLTC883zCfNVMuVa9HgulTp5zPNAYFQcv2AaA6VF/Jg0x+PvSY1I8ek/o5wt6SEuPn79Dqjvyy8HlsB5KwnDxM6OnDRBYcJlynV2gn3T5/uKwQKs43t1CIRReCNZ2yPy75uSUsPetP1Ay28ag9YK3JVXzHNz0f5LHHYMoUj6QhL6fsYPfhw565xtKl0IFUhvOrscNshvGu/T5yhQqw0PXv19L179eC1mRuSeHA/I3kr99GfHgX1FYjoLfZoDU1Jyv3pYRwWzrhheksWRVIYS364s1gugNQ9hamAkOrqxBamM4FW1z7mOzqxNnk4vzl1YMsnsTFYJggLtlcPpi+gCM8xmsu1d9LD67+rXwwfQP7eZA3XKqfwAW8ufbecvtmspV7+K9L9T8lm/+turXcvudYx/m85VL9DXuCePunyeX2jWUF5+PaclRf7+vJu1xabt9NLGE0q1yq/8b+i/nsrG+Fh1ngyNFak38euIVV9Cy37zU+LTfntTq3Jj3KfpxTIYIo4T0+dqkuwOUpb5S7Ug9yeJ1PXaqbQxBjk94tt+8CTvKci1/7sqmq2raFmTPhnnui8PNz/1y4qpgjnJ8rm0+7f55HYWHpUsbGH5wePao/39MCW/rT67o4uC6uwtuGnGzFoV+2kbX9MAWHTlJyLA2dloY5Iw2/02kE5qZhKcomsOQ0Law55FA+XVowtYuGCigfsPng+sozxVT8+NpE9U/428NiwIqPLiYvr/zxEvKJJK3KN55lHSmOZM9ZA4kXk053dtZcGdiZH8Oms/5eTuAgcaV/yGvwy+l+/HxWhrDpbOc8F9/MDGITP3A5ZrMxLfrxx6FP7IfGMH8DKhtM57r2K69Wyj6nFhXl+aCvLnx8oEsX6NKlI9w8s8Lx4tP5pG09Ssa+DPIOpVN0LJ2cM2Ye6Wy8WcjIMF6dDvvz+8HeBFpP08KWQ5DOwVzFz0Qh/hX2Wah5ACffFMiJyL50uvJidr6lMZkb5gvarp0R21qtxhv3ggL3v99bvhyuLTu4c+GFxvx6T1CKlgNjOG9gDDDZ8bu4uBhSU+HwjmDWLHiLwoPH8Uk7hn/mcYJyjhFYeIqg4iyCdTYWe/hcgD+F1O6L4c1gurLvmAq/cpVSM4AZAHVb40yI5iGPAPL9Qxh/6SGGjcwiPj4TPz/N2rU113Wn4wVlAsBTJ0hISCA3N5eEhAT3tL+riAu0L8nEUBwZwfr1G2uu5E1BwHBfGN4B6IACbECB/VXWO/oQBQWp5OX5kJdnpuC0jc93vQ0Fxej8YlRBERQUowqLUIVFmAqLMBUVYi4qwlRcxIWxWZwpLsBqVVitirC8M/zy+2WYbSWYbCWYbSX42Eow6RJ8dAk+tmLM9nIRfkS0LHTUtdnAUmLldEEwZqwYY8rGy9gu/+vaVklyqJqC8bJ0pX8SXFdZfeVKFO+G6+eqIIZEJ/H36/YwYkQGLVsWk5YGCV5I9HLiRBST2cx4FhH7zkl2mkaSduGFdW7v7J/d335rCfQHICAgk4SErZVXPBd0BboGAh0JBC4h4awTAkktOzVAa0pyiijJLKQkq5CiHCu2PCu6oJgCFcB/2myjqEhRUmKiuNhERGpHFhx+DMwmlA8oswmTr8IcEYh/+xZYuoRgjW7jeEo56RfXBpzcJTx8GCdPGkHjggXr6dAh321taw0JCcP5okwwvS8ujqNu+jtQa8Fgm9YdE92xAfn2VymbDfIzrRScKKQwvYj/WLaSk+PLk0+61rzSXsp8rZQaDszSWl9m334YQGv9dFV1uga10+9e/EDVbZYpbx31V2w+zhyG/nmZ9N74YcVKlbCZfdk++s5y+4KyUum6bWEVNezXt3egMCCMvUNuKncs7OQ+ovdUPcJR9t39mdD2JA8on08yIjWRtsnrXOg9nG7djdRel5Tb1zZpLa2ObK+iRvnrZ3Tox8muw8sd77DnJ4LTq07sXvZrn9ZlKKeinFMKlILobYsJyK5kmkaZysePHaNtu3Ycjx1DTuuu5ep32rgA3/zKUzydPTJyJG4c+WHtyh2PWfsJJlsVo3Rn1T88cDzFgaGOuqqkmM7rKx8ZruzPb8qwKWgfX0d9n7zTRG351qW+20w+HBp+fbnjluwTtNm5omJdix8BHcIJjGpFcKdWtO7fDrOf9x9e2/TCKgbdNwaAbaGj6Jf1CwkJCYypbU6kqtqf9S2D/vUHADa0vJTBp5a6pV1RO9qm0TaNrcTGyh9XMmrkaGx+Fmw2nK+8AnT2aWzFVmwlxpxRW4kNXWLFVmO+q0AAACAASURBVGxFa7BZtfHys1DcvpNjepPWYDqVjk/aMcf10MartKw1jn0lwS0p6NDVUVdr8Dt2EN+TRxxzpiq0YT9Z2zSFER3Ij+7hqAsQeGA7fulH7Y354Rfkh28LP/yDjX9D21hoEdO61tkJPOntt+HYjMccDyHy+OMwa1ad2zv7Z3fuXJg+XXM7bxMztgsPv2tPXdMYh6hFtUaOhM5rP6EH+7jj4iTafv0/t2UNSU6GPl3yyKSl84HtI0egfXu3tN9QlFKbtNY1PmzkzZHpDUB3pVQMcASYAtxYXQWfDiGMWVT5nOmzVXwmqSXg+mT+iu/jo4D/c7n+JRX29LC/6mqA/VVXI+yvurrI/qqrK2s8IyEhgdFVBVzTJ9bj2sAtU+tR2Rem31TzaXYjK+wJAepz/TbADfWo37AC24c5y4Xun+aRtzPZWY70zKIwombKpFAmhcnHhG+QLwEtK/lYtKUFOtTns+MI+6uuOtlfdXRhHOC+xUgaQnAw7PVgnuljx4xMHm9xB6wA+odAVpZbryEaRseO8K+1j9ONA/AjkPKw2/KFr10Lo1jtDKR79z7nAuna8FoGVK11CXAXsBTYDXyutXZtcpwQotFqEeWcMx1U7P5gWienOMrWjjFVnyhEMxQcDNmUWXLaA8F0Vw44d3TtKqPS56joaEimzO/QJPctK75mDVzET84dF19c9clNgFfzTGutlwBLvNmH/2/vzuPrrOu8/78+WZvlZGmztWnSje4sBYoUKZAKoiiUQURvRdxQnN846jA64yjjiKPzu11wvJnbfRn3ZZgRRVyQNQhIWYQCpQvd16RNmr1Js37vP06S66TNcvb1/Xw88sh1nXOdKx8eFyd593s+1/crItFVssAL0yXD0Q/T+Ue8ken8FQrTIoGKi09ZAbGzc+qDw9DUBHWBcwcs0KdDqaq+3r9wy7gohuk//xneGximL4/kk+3kp7W5RCSqSuYVM4S/d7uQPoZOhDLB0MzK2r0wXXKOwrRIoHiMTNcHTpgegwVFJD5OG5neu3fqg0PQ2wv7XuxiDZsBcNnZcFkcFgRIoBnDtJkVmtmnzOw7o/tLzezq2JcmIqkoK9voMG90unNf9Ean3Yij5qT3C796ncK0SCCf75SR6SiH6SNHThmZVphOWbEK0y+9BJ2uhApaeV/dH7GvfQ1KS2d+YQoLZmT6+0A/MDa9wyEgyMlCRCQT9WR7NyF2H4zezUkde9spHV1I4QSFVK6qjNq5RdLBaSPTUWzzcG6SNg+F6ZRVXx+bML3ZPyBNJ2X0rr8SPvCBqJw3mQXTM73EOfdWM3sbgHOuz0x3G4jI1F4ovZStx5fSTjmrT+YT4vz3U2p+ci9jY95H8hexNMOXsxY51Wk901Ecme7s9C/uMSFM19dH7fwSX3PmwOHcehj077u9ezHnIr6hdCxMA6yJZBKyFBJMmB4wswJGF1QxsyUQ0iqLIpJhvnru93jwQf/2fQWQz/6onLdj877x7fZStXiInKqoiAkrabqurgiXw/E0+af8Vs90mjCDrCofPYeLKOYE1t3tX/qxIpLpKDMzTAfT5vFp4D6gzsx+CjwETL1yiohkvHKvZZr2KE7o0b/LGxE7WaURMZFTmYEr9tYUt+5u/5rRUdDUBHn0U83ouu1ZWf51qSVlVVUPRLXVY3gYml44RgOPUEQP55wTYYEpYsYw7Zx7AHgT8G7g58Ba51xjbMsSkVQWqzB9rLeYFzmL48xmZN786J1YJI0Ul2TxUe7gb/gabXf+2FvSMUJNTTCfQ94D8+ZBTkJn2JUIVVWdjGqY3rULGvp+zyO8hk5Kqf5c8IvdpbIp3wVmdt4pD41+wEO9mdU7556LXVkikspiFabvnv0+/ov3AfCjG6MTEETSTXEx/DsfBeBDr4XZUcq7zc1gOO7mOtZWHaT+XI1Kp7qqqn7u4Vp2s4SlVy7i6nPPjeh8mzfDOjYBkM1IxnxyMd1b7Muj32cBa4EXAAPOBp4C1se2NBFJVcv7NvNJfkcZHdQ8eQG8uioq5z0UMCg2v043H4pMxud1edDdHb3ztrbCbs7geu7msx+Cf/7n6J1bEqOqqp87uBmAm6rh6qWRnW/zZnjraJgGYN26yE6YIqYM0865DQBm9gvgFufcS6P7ZwIfi095IpKKFrc/y3vw/6V9fMd7GOKdUTnvhDCtLg+RScUqTLe0eNuVmpUyLVRWnhzfPnhwmgOD9PJz/azmZe+BtWsjP2kKCOYGxBVjQRrAObcFyJD7M0UkHLmV3jzTuT3R6fMYGYHDh7392tqonFYk7QSG6Z6e6J1XYTr9VFd7k7NFI0wPvrCVXIYAGKhbAiUlM7wiPQTTSbXNzL4L/AT/9HjvALbFtCoRSWn51V7TdH5fO31ROGfrni5uHfomh6mlzbeQwsKLo3BWkfRTXAyf5N+4hntZ+TedYF+AjRsjPm9gmI5w9jRJElVVE8P0yIh/kpZw9PRAzVFvXryctZkz7hpMmH4P8P8BHxnd/xPwjZhVJCIpr6DWC9MFJzuIxhqIrc/s5Yt8HIDdAyuBrVE4q0j68flgIftYx1NwBP+dg1HQ2gp/z5cpoI+lj9bA6uv8K39IyiooGKasDP6p4+OcN/AcI0v2kvVYY1h9dDt2wBq8MJ11bobMi0cQYdo5dxL4yuiXiMiMiud7YbpoMDptHp1bvR6PzmI1TItMxeeLzSqILS3wN3ydJeyBfwHeconCdBqor4eGjkYu5GnYB+zZE1aY3rZtYpjOmBVbCKJn2sz2mtmeU7/iUZyIpKaSBV6YLhmOTpg+udsL031z1DAtMpXTwnRnZ8TnHByEjg7H3PFZcsmYac/SXV0dUZlrettWxzm84D2QQWE6mDaPwFsxZwE3ALNjU46IpIOS+SWMYGThKKGbkYHIV2AbOeBN5TFcrTAtMpXiYjhIqfdAFEamjx8HH90Ujt0BUVAw8U5HSVmnhek94Y2XHn9uP2X4/+HWX1ROfgZNuRTMCojHA74OO+f+D/CaONQmIikqKyeLTvNm9Og/GvktiNlHvZHp7HqFaZGpxKLNo6WFiaPSNTX+tcsl5dXXR2dk+uDOk9zDRvZTT//KNRn1/8eMI9OnrISYhX+kWv8cFZFpdWWXUz7kb/EYPBZ5mC5o88J0wdLMGfEQCVUs2jxaWqCGgBsZ1eKRNurq4OEIw/TgINx/YAW/5x4Auu4biFZ5KSGYNo8vB2wPAXuBt8SmHBFJFyfyyhmdbpShlt6Iz1fW47V5lK7SyLTIVIqLoTPKbR6nhemamojPKckhGj3Tu3fD0Ojv+/nzwTcnL0rVpYZgwvTNzrkJDTRmtmiqg0VEAJ6cez0P715HB2VUWmQrPDgHlQPeyHTlGoVpkanEYmS6tXWSNg9JC/X1cIB6hskimxE4cgT6+yE/P+hzbN/uba9cGYMik1wwU3P/T5CPiYiMu+/cT/Ahvsqn+ByHcxdEdK6OI73Mxt8yMkAuvsVafk1kKmrzkFDU1sKQ5XGI0fY552D//pDOsS1gKb9MDNNTjkyb2QpgNVBqZm8KeKoE/6weIiJTKvdmx6OnJzeic7W+1MTY6Vpy5lKbHeYSXSIZwOeLTZvHBWrzSEt5eVBdDXubF7GAA/4H9+yBZcuCPkff43/hu3ydzaxhfeHFwHkzviadTNfmsRy4GigDrgl4vBt4fyQ/1MxuAG4HVgKvcs49G8n5RCT5BIbp7u5gOsqmdqw9l0e5mWqOkj2nAjV5iEytuBhaqeCt/IKsslJ+fk/5zC+awaSzeUjaqK/3h+kGHvU/EGLfdOlLj3Mz/wlA0+Z3A9+PboFJbsq/cM65e4B7zOwi59yTUf65W4A3Ad+K8nlFJElEM0wfoJ73810AbrgU3hDR2UTSm88Hg+RxF2+lcAC4MPJztrbCXbyFHSznLZc2U3PGGZGfVJJGXR189+n38Qeu4r3/uojX37g86Nc6B9VN3sqHRRdnzmItY6Zr8/hH59wXgbeb2dtOfd459+Fwf6hzbtvozwj3FCKS5Ja3b+I7fJdy2ml/8UzgM2Gfqzng0+Xq6shrE0lnRUXedm8vDA9DdnZk52xpgUe4GYBL74SaFZGdT5JLXR38kosBOD8PXl8ywwsCHD4Mq4a8lQ996xWmA421kye0BcPMbgFuAaisrKSxsTGR5UgM9fT06PqmkcHdz/I+vgfAw8cGIrq2zzyzGKgHoK9vD42NB6JQoUSL3rvJp6BgPX19/j/x9933GEVF4a1COnZtDx9+NeCf7mznzifp6OiPVqmSQGPXd2BgPuD/tGHTpsM0Nu4M+hzPPVXMB3l5fP+J7i6GMuz3wXRtHveOfv9hOCc2sweByZqqbhttIQmKc+7bwLcBli9f7hoaGsIpR1JAY2Mjur7p47kn++GX/u2i/m4ujODa/uAH3va6dYtpaFgcUW0SXXrvJp/SUugbXStpzZpLqJ3nwlqRrrGxkcsua5hwD+PGjReFMmuaJLGx925rK3z96/7HhodraWgI/s6Uw394kXz8i7S0FC9k/TXXzPCK9DNdm8e9gJvqeefcxulO7Jy7IoK6RCTFFczzmqaLBiObmuviJ+/gfPZzlGoWcSOgqe5FpuPzwb83v43LeYiKxV3wwP1w6aVhnaujw98qMnZeBen0U1fnbR884OBYi//Gl9yZZ2Iaetbrl+5YcA6ZOHHpdG0ed8StChFJO0XzvTDtG+qI6FyvOnQ35+C/D3rb8KUoTItMr7gYSuiiihYYIKLp8Vpa4NU8wSf435ygBn54CbzrXdErVhKu3t9Fxz1s5PIXHoLqXti8Gc45Z8bXFr7ihWl3dub1S8P0bR6Pjm2bWR6wAv9I9Q7nXESLrpvZdcD/BSqB35nZZufc6yI5p4gkF199QJgejmxkurT/6Ph22XLdgSgyk9Pmmo5g4Za2NljKTq7md/7JcR84qTCdZqqr/YPQuYODFNHrf3Dv3qDC9NyjXpguuVRhelJm9kbgm8BuwIBFZvYB59wfwv2hzrlfAb8K9/UikvxK6svGt8voYHhwhOzc0BdbGR6GimEvTM9eqTAtMpPi4lNWQYxgZLqtDSpp8R6ozMQP8tNbVhbMnw+79gZMefjKKzO+rr3NsWrQC9OVr83MMB3MX7YvAxuccw3OucuADcBXYluWiKS67PwcuikGIAtH16Hw/pi3HeihmBMAnCSf/KrSGV4hItEcmW5vV5jOBHV1sJ2AOQ+3b5/xNbsbDzKbdgC6skrJXrwgVuUltWBWUjjmnNsVsL8HOBajekQkjXRll+Mb7vFv72+nfFHZDK84Xdu2o+M3tLTlVDNP89OLzMjni+7IdAWt3gMVFRFUJskqnDD9cvMc/pV7WMNm1qwY5E0Z+vs5mDD9spn9HrgLf8/0DcAzZvYmAOfc3TGsT0RS2Incchg+6N8+1E44Nw52vuK1eHQWVDMvWsWJpLHT2jwi7JleqJHptFdfD4+cGqbd9FMqvrSniHvZyL1s5DNv9S9tnYmCafOYBRwFLgMagBZgNnANcHXMKhORlNc7y7sJse9Ie1jn6Nvnhelen/qlRYJxWptHBCPTavPIDIsWwRHmjbfn0d7un8plGlu3eturVsWwuCQ348i0c+498ShERNLPo0vfz4+fuYZ2yrmuZAXnh3GOgYNemB4oV5gWCUa0R6YVptPfkiUAxnZWcMHY4tfbt0NV1ZSvCQzTK1fGtLykFsxsHouADwELA4+fadEWEZEt59zId5/xb78qzHO4Zi9Mu0qFaZFgRLNnWiPTmWHx6MKyE8L0tm1TLvZzonOII/tHgDyys2Hp0vjUmYyC6Zn+NfA94F5gJLbliEg6Kfe6PGgPr8uD7Jbm8e2suQrTIsEoLo5em0d3az8ldAPgsrOxstBvJJbkV1cHOTmwbShgiHmamxAP/eJxeriSrazi8bK/Ii/v9tgXmaSCCdMnnXP/EfNKRCTtRCNM3++7nkeopZqjrD1/bXQKE0lzPt/YCOPTrN1Qwjd+UT7zi6bS6s3kMVQ6h9ys0OeLl+SXnQ0LF8L2Xf6bEIcLisie5vgTD20ij0HW8AItxRfEpcZkFUyYvtPMPg3cD/SPPeicey5mVYlIWohGmP7D4BW8xBUA/GVDFIoSyQA+H/RSxLNcgA9g6rbXGR3sLOEmfkQFrXz6w4bGpdPXkiXw0K7LWcE2Pv/jpfzV9VPH6fwXnhrf7ll9YTzKS1rBhOmzgJuA1+C1ebjRfRGRKS1r/hMPcTvltHPkkUuBO0M+R7PX5UFNTfRqE0lnxcXednd3+OdxDg51+vgJNwHwb/8QYWGS1BYvhj9SShel7Nk/zYHOMXf/pvHd3EvWxb64JBZMmL4OWOycG4h1MSKSXspyejiPRwAYbA+933loaMInzLrvSSRIPp+33dMT/nlOnsxiYPSvf34+FBREVpckN/+MHn67d09z4MGDzO73j3R0U0ztFRk8lQfBzTP9AuhTHREJ3ay5Xp9HwcnQ+zxaWvwjY+BfdC03N1qViaQ3b2TaMdR5Ao4c8f/rNEQ9Pd6brrx82vU7JA2MzegBsGfP1McNPOa1eDzDBSxfNV13dfoLZmS6GthuZs/g9Uw759y1sStLRNJB0XwvTBcPhB6m2zYf4CnezFGqOcpq4PNRrE4kfY2NTG9nBcubXoFa4JVXQp6/rKvLiwmzZ0exQElKZ5zh/26MMLBlF/z8L3DiBLzvfROO6/jt4+Nt+K+UXchrCuNbZ7IJJkx/OmDbgPXA22JTjoikE1+d96FW8XBHyK/v3n6Ydfgnqt7R3zzD0SIyZmxk+gRF3oNhTI/X3Z3D/88nuJrfMnC4Eu77R3j966NUpSSbZcv8s3osGd7JQ4dWwNvxL9py880TPpbIeezh8e2W1Q3xLzTJzNjm4Zx7FOgE3gj8ALgc+GZsyxKRdFC60BuZLnPtDA+5kF7fu1dLiYuEIzfX3+M8Ya7pMFZB7O7OZSk7OYstnN/5SEQrKUryy8/3j07vZCmdY4v+HDsGu3Z5Bx09yuzDWwAYJIe816xPQKXJZcowbWbLzOxfzGwb8FXgIGDOuQ3Ouf8btwpFJGVlF+bTi/+OpVyG6DxyIqTXDxwKWEp8tqbyEAlFNFZB7O7O0eqHGWb1anBk8ScCVj586KFJtzexjtWvCvj0I0NNNzK9Hf8o9DXOufWjAXo4PmWJSLroyvZGp7v2h9g3HbCU+EiVRqZFQlFcfEqYDnNkWmE6s6xe7f/+EJd7DwYEaDfi2Gv+OxX/yOtYsyae1SWn6cL09UAz8IiZfcfMLsffMy0iErSeXC9MnzgUWpgOXEo8W0uJi4TE54t8SXGNTGeeScP0/fdDXx8A+y6+kcVuF+ewmV+Xv5fa2gQUmWSmDNPOuV85594KrAAagVuBajP7hpldGaf6RCTF9eV7NyH2HgntJsT8Tm9kOr9eYVokFKeNTIcRpk90GXM47j0wZ04UKpNkNhamt3Am+3JGJ57u6oJ77wXg+ecBjBc5h5rz5mm6RIK7AfGEc+6nzrmrgfnAZuCfYl6ZiKSF/gJvZLq/ObSR6eIeL0wXLVaYFgnFaT3T4dw82NZDFv4bh/sLyzTZewZYtgxycgCM7w/d5D3xrW8BsMlb+JDzzotraUkrmEVbxjnn2pxz33LOaSlxEQnKw+d9jDfxSzbwMLtqQrvru7TfC9NlK3QDokgootHmkdvhfZo0VKYWj0yQlwfLl/u3f8Q7cVmjUfHhh+Hpp3nySe/YV786/vUlo5DCdLSY2ZfMbLuZvWhmvzIzrbAokqZaVl3Gr3gTjWygeTD4j4gHBqByxAvT5Ss0Mi0SimjcgJjX5b1mZI7CdKa48EL/930sYvOFHxh/3N1wA01PHxzfv+iieFeWnBISpoEHgDOdc2cDrwCfSFAdIhJj5V6XB+0hdHm07O+lhG4A+skje47+zS0SitPaPHp7Qz5HQY83Mm1VCtOZYt06b/uLlXfA614HgB04wLsH/O0eixdDtcY4gOBWQIw659z9AbubgDcnog4Rib1ww/TR4zm8m/up5igranv4Z93lIhKS4mJ4gNdSQQsfvb2ET3w6L+RzFPW2jW9n11REszxJYoEjzo8+U4g79Ads+zbu/NYsPvUf/mnxLrkkQcUloYSE6VO8F/ivRBchIrExIUy3OYKdYbPpeB4P8loArlwdg8JE0pzPB/3Mop9ZdJ0M/fUjI/Cz/ht4lnOooJX/+aDmQMsUK1dCaam/M6ipCba8bJx11ip+FnDz4VVXJa6+ZBOzMG1mDwKT3TF0m3PuntFjbgOGgJ9Oc55bgFsAKisraWxsjH6xkhR6enp0fdOQe2QXe/kcZXTw4gPraGwMbjKgRx+twT8zJzjXTGPj9hhWKZHQezc5NTfXAksB2L79MI2NO0N6fXd3Di2sp4UqCguHeLz/cdB1TivTvXfPP38lDz/s7+O48849XHVVM888cxFgZGU5CgqeoLFxKH7FJrGYhWnn3BXTPW9m7wKuBi53zrlpzvNt4NsAy5cvdw0NDdEsU5JIY2Mjur7pZ9uOLBayH4CS4W7WBHmN//xnb3vNmhoaGjSbR7LSezc57d3rbZeU1NLQENrI8u7d3nZlZY6ucRqa7r3b1OSfwAPg+ecXs3DhYsbS2iWXGBs3hjY7UzpLSJuHmb0e+DhwmXMu9DsiRCRlFNZ6fR5FA8E3TR9t9lpCdJOLSOh8Pv/3Og4w+3AnPNHln8ssyPsPAu9xCGzXksxw1VX+afIGBuC55/xfY265JXF1JaNEzebxVcAHPGBmm83smwmqQ0RizFfv/RX2DQUfpq/6w4c5RiUvchZrD9wdi9JE0lpxsf/7NlbylYfOhvXr4cSJoF/f1ga5DAAwe3YsKpRkVlYGN954+uOVlXD99fGvJ5klajaPMxLxc0Uk/koWeGG6zLUzMgJZQfwzflZHE5W0UkkrW3wDMaxQJD2NjUx3UUIRox8Cd3V5KXsG7e2wl0WU0smJTZXQ/CTUqN0qk3zmM3D33ROnKL/jDsjPT1xNyShRI9MikiFySgoZwL8E8Sz66WzuC+p1vp7m8e2iJfoDLhKqscw8YRXEEBZuaTvuqKCVYk5Q3bsPSkpmfI2kl7o6eOQReO1r4YIL4Hvfg3e+M9FVJZ9kmBpPRNKZGZ1Z5VSOHAOg60AH5fMKZnxZ2UDAUuIrFaZFQhU4Mj0uhCXFTzR3kz/a5jGQW0heYWE0y5MUce65cP/9Mx+XyTQyLSIx15XtrV7Yve/4jMf39kLViDcyXbpUdyCKhCrSMD14pGV8u69Iqx+KTEVhWkRirivfu3upZ//MYfrY3hP46AH8S4lnzdZS4iKhirTNY+SoF6b7SxWmRaaiMC0iMdc7ywvDJw+1znh82zavxaMtryboqbxExDNrFmRnhz8yTav3Xh2erTAtMhWFaRGJub4iL0wPNs8cprt3ei0eXYXqlxYJh5l/dDrckensdm9kmjkVUaxMJL0oTItIzD14/vs5l+eo4wB/WvTuGY/v2+uF6b4S9UuLhMvnO2VkOoQwnd/phensGo1Mi0xFs3mISMwNLJjPZpYBcLRj5uOHDnttHoNzNDItEq7iYugg4J6DEMJ0QY8XpnPnKUyLTEVhWkRirrR0cHy7deYuD2j2RqapVpgWCZfP5w/TPRSRW1VOflFR0K8t7vPC9Kw6hWmRqajNQ0RiLtQw/cOaj3MGO7mYxzl+9btiWJlIevP54Pu8Bx89PPbTg/C5zwX1uv5+KBv23qwK0yJT08i0iMRcaekgxgjltJPbdAKon/b4A62F7OYMdnMGvjXxqVEkHfmnx/PPhtPdHfzr2tvhen5JBa0sLWul8bIFMalPJB0oTItIzM3r28sg68hmhP17FgF7pj0+sMujRl0eImEbW7gFoKcn+Ne1t8MA+RyhlqLKWtBU7yJTUpuHiMRc/txCshkBoHyklZGRqY91Do569x9Srck8RMI2tnALhDYy3dbmbc+ePfVxIqIwLSLxUFrIENkAlNBN+9GBKQ/t7HDM699DAb0UF08MAyISGp8PjBEu4U9Ub7oHfvKToF7X3u5tl5fHqDiRNKEwLSKxl5VFR/ac8d32XVMvKX5sZyd7WEIvRezumxuP6kTS1libRyMNXP/jv4KbboLBwelfBHQ297GSrVTQwpyy4RhXKZLaFKZFJC668rwV1Lr3Tj2lR/v2gDmmcwtjWpNIuisuBkdWyKsgZm17ma2spoUq/veDF8SwQpHUpzAtInHRW+CNTPcemDpM9+zy7j7sLtLdhyKRGBuZbiegVyOwh2MKQ03eHNMDvjnTHCkiCtMiEhf9xd7I9MnDU4fpvn3eyPTJUoVpkUiMhekJqyB2zLwMqbUcG98eLNMc0yLTUZgWkbgYLPPC9GDT1D3TgweaxrdHqjSVh0gkSke7O0Idmc467o1Mj1RURbsskbSiMC0i8VHpjW655qNTHpbdfHh8O6tufkxLEkl3JSX+76GOTOd1emHaqjQyLTIdhWkRiYucWq9lI7ulecrjZh0/NL6dv7g2pjWJpLuxMB3qyPSsbi9M58xTmBaZjsK0iMTFrIX+MD1ALgMnpp5nuqTbG5kuWakwLRKJScN0ECPTRb1ez/Ss+QrTItNJSJg2s8+a2YtmttnM7jezeYmoQ0TiJ/ev3sgcWpnFSf6u+HuTHjMyApX93sj0nLMVpkUiMdYzPaHNI4iR6ZJ+b2S6cKF6pkWmk6iR6S855852zq0Bfgv8S4LqEJE4qVlcSBtzcGTR1ORfNvxUx446avFGpmedoZ5pkUhMOjXeDCPTzkH5kBemfYs1Mi0ynYSEaedcV8BuETDJn1URSSclJTBrln+7txe6u08/pumVrvHFJXqyfF4SEJGwZGdDUREcYj4vcDZDF18KS5ZM+5rubqjEC9N5tQrTItPJSdQPNrN/A94JdAIbpjnuFuAWgMrKShobG+NSn8RfT0+Prm+a6unp4dFHGykvv5Cm5P3uOQAAEpFJREFUpgIAfv3rp6iv75tw3OOPz+FTHCWPfq48Zycf1f8PKUHv3eRWUHAR957YyL1s5K6PPEllZT9Mc72OHcrmYkrIpx/D8cTmzWAWv4IlbvTejY6YhWkzexCYbMWF25xz9zjnbgNuM7NPAH8LfHqy8zjnvg18G2D58uWuoaEhRhVLojU2NqLrm57Gru35tc00Nb1MDc0sLns16xtmTzju5Zf93wfIZ+7aM9H/DqlB793kVlEBraPrJK1efRGrVk1//PPPw3wOA46zFx7hhQ26dyFd6b0bHTEL0865K4I89GfA75giTItI+vjCvrewiscAePipB2Hj5ROeP+y1S1Orv98iUTE2owdAZ+fMx7e1jW0ZWWWFsShJJK0kajaPpQG7G4HtiahDROLrZNnc8e3+fU2nPX/Im8iD+br3UCQqAsN0V9fUx43xwjT4fEPRL0gkzSSqZ/rzZrYcGAH2A3+doDpEJI6GqubCLv/28OHTF24p2fJnLqeXQ8ynrnIhMCuu9Ymko7Hp8d7Gz6j++TF4rBVuvx1yJo8ACtMioUlImHbOXZ+InysiiZU117uNIuvY6SPTb972WRq4D4BDB38NXBuv0kTS1tjI9J18hMofjjZPf/jDUDX5/NFZ217mDeyjhUrm5eUBcyc9TkT8tAKiiMRN3gLvj3J+28SR6cFBqD65b3y/cu2CeJUlktbGwnQrFd6DY3ckTmLJpp/wO67maS7kjQd/FOPqRFKfwrSIxE3JMm9kurBr4sj0wQOOhewb389fsSheZYmktVDDdF7b0fHtwfKyKY8TET+FaRGJm6rzvLsKq/oOMDLiPXfo2WYKOAlAV0651+gpIhEZeysdZ4734DRhuqDL+9RoqGL2lMeJiJ/CtIjETeFKr3WjjgMcOTg8vt/+/N7x7eMlGpUWiZZQR6ZLTnifGg1XlU95nIj4KUyLSPwUF9Oe4/+Dnscgh5/1/mif3LZvfLu3amGcCxNJX6GG6bKT3sh01jx9QiQyE4VpEYmrwFHn9uf3jW+P7PZGpt0CjUyLRMtYm0dQYXp4mNlDx8Z3c2pLJj9ORMYlap5pEclQ3TVL2d12nH0s5HBz9vjjsw7vHt8uPkthWiRaQhqZbmkhm5HR4+dQPNtiXJ1I6tPItIjE1dMf+glnsJsreIg/DV4E+KfFm9u5bfyYqkuWJ6o8kbQTSpgeOOC1eDRTQ0HB8KTHiYhHYVpE4mrZcm+ka+tW//fduxwrnBemC9euindZImkrlDDds9O7j+F47lxMA9MiM1Kbh4jE1VlnedtbtsDwMGzffJKn2chKtrEwv4nKuVpxTSRaxnqmD1PL77Ov4Q03zYHlk3/607fPG5nuKKhBtx+KzExhWkTiqqIC5s6FpiY4eRJ274bNOwr4DP6V1j7yAcf/0XCYSNT4fGAGh1wdbxz+DQPfhtzcyY/tsjJ2cwk1NHPct0hhWiQICtMiEndvqd/EUNOznMkWXnnkU/z5z7Xjz13wKgVpkWjKyoKyMmhv9+93dEBl5eTH7jzzOq7lOgDeuAY+RmN8ihRJYQrTIhJ3HzryTyzhUQB+8MDVPPmkF6bXr09UVSLpa/ZsL0y3t08dptvavO1yrdciEhTdgCgicTeyymuc3vebF+np8W/X1kJ9fYKKEkljgcE4MDCfaixwgz+Ai8jMNDItInFX/bo18Ef/9rWD/82lPMjjrIelGzDbkNjiRNLQWJh+M/9NzeefgPwmuPVWWLduwnGBk3woTIsER2FaROKu5JrL4O/92+eyGYDX8AhHjz0OKEyLRNtYmL6a37LwHv/Nvlx55cQw7Rzrf/9J/pZ5HKCeyoqN8S9UJAUpTItI/C1ZQlfVEkqO7Z7wcPlbrkhQQSLpbSxMNxEw7WRT08SDjh/nqs2f5yqgCx9/rOqKW30iqUw90yISf2b4/vodEx5yWVnkvevtCSpIJL2NtWwcYZ734JEjEw86eHB88wD1VFQgIkFQmBaRhLBb/w4WLvT2/+EfJuyLSPQENTIdEKYPUjfljB8iMpHaPEQkMcrKYNMm+MEP/CH6hhsSXZFI2goqTB84ML55kDrOrZhy1XERCaAwLSKJU10NH/94oqsQSXtjYXq6No+R/QfHP64+QD1z5sSnNpFUl9A2DzP7mJk5M1NnloiISIyM9UxPGJluboaRkfHdgT1em8fxwroplxwXkYkSFqbNrA54LXBgpmNFREQkfGMj0ycpoDN7dGdwEI4eHT9mZJ/357inrC6e5YmktESOTH8F+EfAJbAGERGRtBe4AuL+rEXezt6945s5+7ypKnurA44RkWklJEyb2UbgsHPuhUT8fBERkUwSGKZ3jSz2dvbs8X/v7iavrRmAfvIYrq2PY3UiqS1mNyCa2YNAzSRP3QZ8ErgyyPPcAtwCUFlZSWNjY7RKlCTT09Oj65umdG3Tm65v8nMOsrIuY2TE+K/hN7PmxkIGamvoyMnhZGMjxTt3snb02D0spn/oGI2NO3Rt05yub3SYc/HtsjCzs4CHgN7Rh+YDR4BXOeeap3vt8uXL3Y4dO2JcoSRKY2MjDQ0NiS5DYkDXNr3p+qaGykpvqrvmZv9kOuP27+ex936fAw/vpIm5tPzjHXzhC7q26U7Xd3pm9hfn3NqZjov71HjOuZeAqrF9M9sHrHXOaTZLERGRGCkv98J0W9spYXrBAn51zu185WH/7hc1x5ZI0LQCooiISAYI7Jtubz/9+cAFWrT6oUjwEr5oi3NuYaJrEBERSXdjc00DHD8+ujE87G+ozsmhpcV7vkIj0yJB08i0iIhIBqiq8rarf3wHvPrV4PPBo48CKEyLhElhWkREJAMEhum8fa/Ak09CXx989auwbBm3bb2RG/kJADWTzcUlIpNKeJuHiIiIxF5gmN7pO481Yzu//jUA17GTbrL5Ke9QmBYJgUamRUREMkBgmH4p7/xJj3mKCykrg1mz4lSUSBpQmBYREckAgWH6meHzJp2y4/e8gblz41iUSBpQmBYREckAgdm5uSUb3v72Cc8/yqXsY5FaPERCpDAtIiKSAQJHpo8dA267DZYtA6C/sIwP8x+Abj4UCZXCtIiISAYIHJk+dgxcRSU8/zw88ADf/OguXuQcALV5iIRIYVpERCQDFBT4p5UGGBqCjg6gsBCuuIK9XXPGj9PItEhoFKZFREQyRGCrx9Gj3nZTk7etMC0SGoVpERGRDFFb620fOuRtHzzobdfVxa8ekXSgMC0iIpIh5s/3tgPD9P793vaCBfGrRyQdKEyLiIhkiMBR57HR6IEBr83DbOLotYjMTGFaREQkQ0wWpg8eBOf82/PmQV5e/OsSSWUK0yIiIhkisM1jLEwfOOA9phYPkdApTIuIiGSIwJHpsZ5p9UuLREZhWkREJENMNjKtMC0SGYVpERGRDFFZCbNm+bc7O6G9HV55xXt+0aLE1CWSyhSmRUREMoQZLF/u7b/8MmzZ4u2feWb8axJJdQrTIiIiGWT1am/7hRdg+/bJnxOR4ChMi4iIZJDA0ed77vHPMw3+furS0sTUJJLKFKZFREQySODo8wMPeNtq8RAJT0LCtJndbmaHzWzz6NcbElGHiIhIppmqleOss+Jbh0i6yEngz/6Kc+6OBP58ERGRjLNoEVRVwbFjEx+/7LLE1COS6tTmISIikkGysuB1r5v4WEEBXHppYuoRSXXmnIv/DzW7HXg30AU8C3zUOdc+xbG3ALcAVFZWnn/XXXfFqUqJt56eHoqLixNdhsSArm160/VNPVu3lvDBD543vr9x42FuvXXnacfp2qY3Xd/pbdiw4S/OubUzHRezMG1mDwI1kzx1G7AJaAUc8FlgrnPuvTOdc/ny5W7Hjh1RrVOSR2NjIw0NDYkuQ2JA1za96fqmpu98B770JTj/fPjWt6Ck5PRjdG3Tm67v9MwsqDAds55p59wVwRxnZt8BfhurOkREROR073+//0tEIpOo2TzmBuxeB2yZ6lgRERERkWSVqNk8vmhma/C3eewDPpCgOkREREREwpaQMO2cuykRP1dEREREJJo0NZ6IiIiISJgUpkVEREREwqQwLSIiIiISJoVpEREREZEwKUyLiIiIiIRJYVpEREREJEwK0yIiIiIiYTLnXKJrCJqZdQM7El2HxEwF0JroIiQmdG3Tm65v+tK1TW+6vtNb4JyrnOmgRK2AGK4dzrm1iS5CYsPMntX1TU+6tulN1zd96dqmN13f6FCbh4iIiIhImBSmRURERETClGph+tuJLkBiStc3fenapjdd3/Sla5vedH2jIKVuQBQRERERSSapNjItIiIiIpI0kjJMm9nrzWyHme0ys3+a5Pl8M/uv0eefMrOF8a9SwhHEtf17M9tqZi+a2UNmtiARdUp4Zrq+Ace92cycmeku8hQRzLU1s7eMvn9fNrOfxbtGCV8Qv5vrzewRM3t+9PfzGxJRp4TOzP7TzI6Z2ZYpnjcz+4/Ra/+imZ0X7xpTXdKFaTPLBr4GXAWsAt5mZqtOOexmoN05dwbwFeAL8a1SwhHktX0eWOucOxv4H+CL8a1SwhXk9cXMfMCHgafiW6GEK5hra2ZLgU8AFzvnVgN/F/dCJSxBvnf/GbjLOXcu8L+Ar8e3SonAD4DXT/P8VcDS0a9bgG/Eoaa0knRhGngVsMs5t8c5NwD8Arj2lGOuBX44uv0/wOVmZnGsUcIz47V1zj3inOsd3d0EzI9zjRK+YN67AJ/F/4+kk/EsTiISzLV9P/A151w7gHPuWJxrlPAFc30dUDK6XQociWN9EgHn3J+AtmkOuRb4kfPbBJSZ2dz4VJcekjFM1wIHA/YPjT426THOuSGgE5gTl+okEsFc20A3A3+IaUUSTTNeXzM7F6hzzv02noVJxIJ57y4DlpnZE2a2ycymGwmT5BLM9b0deIeZHQJ+D3woPqVJHIT6t1lOkYwrIE42wnzqlCPBHCPJJ+jrZmbvANYCl8W0Iommaa+vmWXhb8t6d7wKkqgJ5r2bg/9j4gb8nyg9ZmZnOuc6YlybRC6Y6/s24AfOuS+b2UXAj0ev70jsy5MYU6aKUDKOTB8C6gL253P6x0njx5hZDv6PnKb7CEOSQzDXFjO7ArgN2Oic649TbRK5ma6vDzgTaDSzfcA64De6CTElBPt7+R7n3KBzbi+wA3+4luQXzPW9GbgLwDn3JDALqIhLdRJrQf1tlqklY5h+BlhqZovMLA//jQ6/OeWY3wDvGt1+M/Cw04TZqWDGazvaBvAt/EFaPZepZdrr65zrdM5VOOcWOucW4u+J3+icezYx5UoIgvm9/GtgA4CZVeBv+9gT1yolXMFc3wPA5QBmthJ/mG6Ja5USK78B3jk6q8c6oNM515ToolJJ0rV5OOeGzOxvgT8C2cB/OudeNrN/BZ51zv0G+B7+j5h24R+R/l+Jq1iCFeS1/RJQDPz36D2lB5xzGxNWtAQtyOsrKSjIa/tH4Eoz2woMA//gnDueuKolWEFe348C3zGzW/G3ALxbg1ipwcx+jr/9qmK05/3TQC6Ac+6b+Hvg3wDsAnqB9ySm0tSlFRBFRERERMKUjG0eIiIiIiIpQWFaRERERCRMCtMiIiIiImFSmBYRERERCZPCtIiIiIhImBSmRURERETCpDAtIpJEzGyOmW0e/Wo2s8MB+3+O0c8818y+O83zlWZ2Xyx+tohIqku6RVtERDLZ6EInawDM7Hagxzl3R4x/7CeBz01TU4uZNZnZxc65J2Jci4hIStHItIhIijCzntHvDWb2qJndZWavmNnnzexGM3vazF4ysyWjx1Wa2S/N7JnRr4snOacPONs598Lo/mUBI+HPjz4P/uXCb4zTf6qISMpQmBYRSU3nAB8BzgJuApY5514FfBf40OgxdwJfcc5dAFw/+typ1gJbAvY/BnzQObcGuAToG3382dF9EREJoDYPEZHU9IxzrgnAzHYD948+/hKwYXT7CmCVmY29psTMfM657oDzzAVaAvafAP7dzH4K3O2cOzT6+DFgXvT/M0REUpvCtIhIauoP2B4J2B/B+92eBVzknOtjan3ArLEd59znzex3wBuATWZ2hXNu++gx051HRCQjqc1DRCR93Q/87diOma2Z5JhtwBkBxyxxzr3knPsC/taOFaNPLWNiO4iIiKAwLSKSzj4MrDWzF81sK/DXpx4wOupcGnCj4d+Z2RYzewH/SPQfRh/fAPwuHkWLiKQSc84lugYREUkgM7sV6HbOTTfX9J+Aa51z7fGrTEQk+WlkWkREvsHEHuwJzKwS+HcFaRGR02lkWkREREQkTBqZFhEREREJk8K0iIiIiEiYFKZFRERERMKkMC0iIiIiEiaFaRERERGRMP0/lLP8juwI55kAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Nlam = 12\n",
    "dx = vs0 / (Nlam * 2 * f0)\n",
    "dz   = dx    # grid point distance in z-direction (m)\n",
    "\n",
    "op = 2       # order of spatial FD operator\n",
    "\n",
    "# calculate dt according to CFL criterion\n",
    "dt = dx / (np.sqrt(2) * vs0) \n",
    "\n",
    "vy = FD_2D_SH_JIT(dt,dx,dz,f0,xsrc,zsrc,op)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While, the direct SH wave is modelled very accurately, the fit of the boundary reflection waveforms can be improved. Let's try more grid points per minimum wavelength ($N_\\lambda = 16$)..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "nx =  1103\n",
      "nz =  1103\n",
      "nt =  2027\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtMAAAFNCAYAAADCcOOfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xd8VFX6+PHPmUx6QhIg1FBCS4CEIlVDFVTsqLgWdBdWxbW7VnTVVVTEsrvWXVdWZH9fFVkL2FBshN57DT100kkvM3N+f9xkbkLakMxkQvK8X695cW455z7kpjxz5txzlNYaIYQQQgghxLmzeDsAIYQQQgghzleSTAshhBBCCFFHkkwLIYQQQghRR5JMCyGEEEIIUUeSTAshhBBCCFFHkkwLIYQQQghRR5JMCyFEI6GUylVKdfN2HEIIIVwnybQQQriZUmqEUmqVUuqMUipDKbVSKTWktnpa6xCt9cGGiFEIIYR7WL0dgBBCNCVKqRbAd8A9wP8AP2AkUOTNuNxBKeWjtbZ7Ow4hhGhMpGdaCCHcqxeA1nqe1tqutS7QWv+ktd4GoJT6o1Jqt1IqUym1WCnVpayiUkorpXqUlq9QSu1SSuUopY4rpR4r3T9GKXVMKfWEUipFKXVSKTWx9Py9pT3hT5dr018p9aZS6kTp602llH+540+UtnFCKXXnWTHMVUr9Sym1SCmVB4xVSl2plNqslMpWSh1VSj1frq2upfWnlh7LVEr9SSk1RCm1TSmVpZR617NffiGEaFiSTAshhHvtBexKqf8qpS5XSkWUHVBKTQSeBq4HIoHlwLxq2vkQuFtrHQrEAb+VO9YOCAA6As8Bs4HbgEEYveDPlRt7/RdgODAA6A8MBZ4pjWcC8AgwHugBjK4ijluBl4FQYAWQB/weCAeuBO4p/X+VNwzoCdwEvFkaw3igL/A7pVRV1xFCiPPSeZdMK6XmlPbG7HBDW2OVUlvKvQqr+KMghBAu01pnAyMAjZHkpiqlvlFKtQXuBl7RWu/WWtuAmcCA8r3T5ZQAfZRSLbTWmVrrTWcde1lrXQJ8BrQG3tJa52itdwI7gX6l504GZmitU7TWqcALwO2lx34HfKS13qm1zi89dravtdYrtdYOrXWh1jpRa729dHsbxpuBs5PjF0vP/Qkj+Z5Xev3jGG8gBrr21RRCiMbvvEumgbnABHc0pLVeorUeoLUeAFwM5AM/uaNtIUTzVZosT9FaR2H0KnfA6KHtArxVOtwhC8gAFEYP89luAK4AkpVSS5VSF5Y7ll5u7HJB6b+nyx0vAEJKyx2A5HLHkkv3lR07Wu5Y+XKV+5RSw5RSS5RSqUqpM8CfMJL58s6OpbrYhBDivHfeJdNa62UYf4CclFLdlVI/KqU2KqWWK6Vi69D0JOCH0t4ZIYRwC631HoxOgDiMxPRurXV4uVeg1npVFfXWa62vBdoACzEeZqyLExhJfJnOpfsATgJR5Y51quq/cNb2p8A3QCetdRjwPsYbAiGEaJbOu2S6Gh8AD2itBwGPAf+sQxs3U/3YRSGEcIlSKlYp9ahSKqp0uxNwC7AGI/F8SinVt/RYmFLqxira8FNKTVZKhZUO5cgG6jqLxjzgGaVUpFKqNcYY649Lj/0PmKqU6q2UCio9VptQIENrXaiUGooxploIIZqt835qPKVUCHAR8LlSzs4R/9Jj1wMzqqh2XGt9Wbk22gPxwGLPRiuEaAZyMB7Ae0QpFQ5kYUyV97jWOrv0d9ZnpeOkzwA/A59X0c7twLtKKR8gCeMBw7p4CWgBbCvd/rx0H1rrH5RSbwNLAAfwYul1a5rG717gb6WzcizFSMjD6xibEEKc95TWZ3+C1/gppboC32mt40rndE3SWrevR3sPAX211tPcFKIQQpx3lFK9gR2Af+kDkkIIIWpx3g/zKH1y/lDZR6XK0P8cm7kFGeIhhGiGlFLXlQ4riQBeBb6VRFoIIVx33iXTSql5wGogpnThgjswpn66Qym1FWNKqGvPob2uGA/dLHV/tEII0ejdDaQCBzDGZd/j3XCEEOL8cl4O8xBCCCGEEKIxOO96poUQQgghhGgsJJkWQgghhBCijs6rqfHCw8N1jx49vB2G8JC8vDyCg4O9HYbwALm3TZvc36ZL7m3TJve3Zhs3bkzTWkfWdt55lUy3bduWDRs2eDsM4SGJiYmMGTPG22EID5B727TJ/W265N42bXJ/a6aUSnblPBnmIYQQQgghRB1JMi2EEEIIIUQdSTIthBBCCCFEHZ1XY6aFEEIIIdytpKSEY8eOUVhY6O1QGlRYWBi7d+/2dhheFxAQQFRUFL6+vnWqL8m0EEIIIZq1Y8eOERoaSteuXVFKeTucBpOTk0NoaKi3w/AqrTXp6ekcO3aM6OjoOrUhwzyEEEII0awVFhbSqlWrZpVIC4NSilatWtXrUwlJpoUQQgjR7Eki3XzV995LMi2EEEII4WU+Pj4MGDDA+Tp8+DCJiYmEhYUxcOBAYmJiGDVqFN999129r3X48GHi4uJqPW/mzJkVti+66KJ6X7spkjHTQgghxHlCazh1Ctq2BYt0hzUpgYGBbNmypcK+w4cPM3LkSGcCvWXLFiZOnEhgYCDjxo3zeEwzZ87k6aefdm6vWrXK49c8H8mPohBCCHEeSE2F+3r9TH6H7nzZ8i4O7cz3dkiigQ0YMIDnnnuOd999t9KxpUuXOnu1Bw4cSE5ODlprHn/8ceLi4oiPj2f+/PmV6s2dO5f777/fuX3VVVeRmJjI9OnTKSgoYMCAAUyePBmAkJAQgGrbLVtRcdKkScTGxjJ58mS01p74UjQqXuuZVkoFAMsA/9I4vtBa/9Vb8QghhBCN2QMPwIL9o3gcBzee+Q/fXOZL9LF/ejss4SZliStAdHQ0CxYsqPK8Cy64gNdff73S/jfeeIP33nuPhIQEcnNzCQgI4KuvvmLLli1s3bqVtLQ0hgwZwqhRo1yKZ9asWbz77ruVesuBGtvdvHkzO3fupEOHDiQkJLBy5UpGjBjh6pfhvOTNnuki4GKtdX9gADBBKTXci/EIIYQQjdK+fTB/Pjiw8CLPAnDp8Tkc3JTl5ciaHqU896pJ2TCPLVu2VJtIA9X29CYkJPDII4/w9ttvk5WVhdVqZcWKFdxyyy34+PjQtm1bRo8ezfr16+vz5QGosd2hQ4cSFRWFxWJxjv1u6ryWTGtDbummb+mr6X8WIIQQQpyjjz82/rVh5UHeBiCAInbO+NKLUQlv2Lx5M7179660f/r06fznP/+hoKCA4cOHs2fPHpeGWFitVhwOh3PblSniamrX39/fWfbx8cFms9Xa3vnOq2OmlVI+SqktQArws9Z6rTfjEUIIIRqjb78tKyk2973duT9k9c9eiUd4x7Zt23jxxRe57777Kh07cOAA8fHxPPnkkwwePJg9e/YwatQo5s+fj91uJzU1lWXLljF06NAK9bp27cqWLVtwOBwcPXqUdevWOY/5+vpSUlJS6VqutNuceHU2D621HRiglAoHFiil4rTWO8qfo5SaBkwDiIyMJDExseEDFQ0iNzdX7m8TJfe2aZP761n5aTbyN0cBMVgs0HJSG9hpHOuRsopff12Kj49nPthtLvc2LCyMnJwcALKzPXed0kvUcLziCfn5+Sxfvpz+/fuTn59PZGQkr776KkOHDq107muvvcby5cvx8fEhJiaGESNG4Ofnx9KlS4mPj0cpxQsvvEBwcDBpaWk4HA7sdjv9+vUjKiqKvn370qdPH+e1cnJymDJlCnFxcfTv358PP/zQGeP48eOrbDc/Px+bzeaMrbi4mMLCwkqxNkaFhYV1/l5XjeUpS6XUX4E8rfUb1Z0TExOjk5KSGjAq0ZDKngIWTY/c26ZN7q9nrfnrDwyfcQXptOSHtlO57dgs8vzCCdZ5AOxYdIS4yzt55NrN5d7u3r27yqETTZ0sJ26q6ntAKbVRaz24trpeG+ahlIos7ZFGKRUIjAf2eCseIYQQojHKXmzM7duKDLpE2cBq5WCrIc7jqT9Xnm1BCNFwvDlmuj2wRCm1DViPMWa6/sv6CCGEEE1Iy93mQhnB4y4EIK+ruXpd4aZdDR6TEMLktTHTWuttwEBvXV8IIYRo7OzFdmKzzWfzO91kLOfsE98bNhj7/A9IMi2EN8kKiEIIIUQjdeTXfYRgjI0+ZWlP5AXG2OiIi/o4z2mdKsm0EN7k1dk8hBBCCFG9Uz9vJ7q0fKxlP9qVljuM78MahnGIaHaVxNHXDj4+3opSiOZNkmkhhBCikSpYt91ZzusW7ywHdW3DxLZrOH0acMCdx6BLFy8EKISQYR5CCCFEYxW4f5uz7D+4X4Vj0dFm+dChhopIeNKCBQtQSrFnT/0mN5syZQpffPFFjefMnDmzwvZFF11Up2s9//zzvPFGtbMauyQxMZGrrrqqxnOysrL45z//6dw+ceIEkyZNqtd13UWSaSGEEKKRap9m9ky3GRdf4Zgk003PvHnzGDFiBJ999pnHr3V2Mr1q1apqzmwczk6mO3ToUOsbhoYiybQQQgjRCGUezaWr/SAANnzofFnFBSUkmW5acnNzWblyJR9++GGFZLps4ZxJkyYRGxvL5MmTKVtwb8aMGQwZMoS4uDimTZvG2Qvx/frrr1x33XXO7Z9//pnrr7+e6dOnU1BQQEJCApMnTwYgJCTEed5rr71GfHw8/fv3Z/r06QDMnj2bIUOG0L9/f2644Qby8/Nr/P98/vnnztUTR40aBRirDE6dOpX4+HgGDhzIkiVLKtU7u6c7Li6Ow4cPM336dA4cOMCAAQN4/PHHOXz4MHFxcTW2O3fuXK6//nomTJhAz549eeKJJ2q5C3UjybQQQgjRCCVvSGUTAynEn2T/XliD/SscH8hmZvIUHzGFnr/8y0tRCndZuHAhEyZMoFevXrRs2ZJNmzY5j23evJk333yTXbt2cfDgQVauXAnA/fffz/r169mxYwcFBQV8913F5Touvvhidu/eTWpqKgAfffQRU6dOZdasWQQGBrJy5Uo++eSTCnV++OEHFi5cyNq1a9m6daszAb3++utZv349W7dupXfv3s7lxaszY8YMFi9ezNatW/nmm28AeO+99wDYvn078+bN4w9/+AOFhYUufX1mzZpF9+7d2bJlC6+//nqFYzW1u2XLFubPn8/27duZP38+R48edel650KSaSGEEKIR2p4bzSA2EUIus8b/Wul4dNEenmIWU/gvXQ/+5oUIm7DnnwelXHtNm1a5/rRpFc95/vlaLzlv3jxuvvlmAG6++WbmzZvnPDZ06FCioqKwWCwMGDCAw4cPA7BkyRKGDRtGfHw8v/32Gzt37qzQplKK22+/nY8//pisrCxWr17N5ZdfXmMcv/zyC1OnTiUoKAiAli1bArBjxw5GjhxJfHw8n3zySaVrnS0hIYEpU6Ywe/Zs7HY7ACtWrOD2228HIDY2li5durB3795avza1qandcePGERYWRkBAAH369CE5Obne1zubzOYhhBBCNEL79hn/2rES2a99peMhvTqY5ZwTDRWW8ID09HR+++03duzYgVIKu92OUorXXnsNAH9/81MJHx8fbDYbhYWF3HvvvWzYsIFOnTrx/PPPV9nLO3XqVK6++moCAgK48cYbsVprTv201iilKu2fMmUKCxcupH///sydO5fExMQa23n//fdZu3Yt33//PQMGDGDLli2VhqFUxWq14nA4nNuu9FzX1G5VXzt3k55pIYQQohEq32HXq1fl4+G9zQS7ZeHJBohIeMoXX3zB73//e5KTkzl8+DBHjx4lOjqaFStWVFunLMls3bo1ubm51T6M16FDBzp06MBLL73ElClTnPt9fX0pKSmpdP6ll17KnDlznGOiMzIyAMjJyaF9+/aUlJRUGhpSlQMHDjBs2DBmzJhB69atOXr0KKNGjXLW3bt3L0eOHCEmJqZCva5duzqHuGzatIlDpQ8EhIaGkpOTU+W1XGnXkySZFkIIIRqhsp5pqDqZbhVnJtNtHScoLqq910+46PnnQWvXXh98ULn+Bx9UPKeWYR7z5s2r8KAgwA033MCnn35abZ3w8HDuuusu4uPjmThxIkOGDKn23MmTJ9OpUyf69DFXzpw2bRoXXnih8wHEMhMmTOCaa65h8ODBDBgwwPkw4IsvvsiwYcO45JJLiI2NrfH/A/D4448THx9PXFwco0aNon///tx7773Y7Xbi4+O56aabmDt3boWe47L/d0ZGBgMGDOBf//oXvUq/+Vu1akVCQgJxcXE8/vjjFeq40q4nKVe63BuLmJgYnZSU5O0whIeUPbEsmh65t02b3F/303YHfw16nZ3FPdhHT349HU9km8ofvedaQgnRuQAc3ZZBp/gIt8bRXO7t7t276d27d+0nnqfuv/9+Bg4cyB133FFhf05ODqGhoV6KqnGp6ntAKbVRaz24troyZloIIYRoZFI2HWNGsTElWZpqTavI1CrPS/PrQEiRMR4kc8cJtyfT4vw3aNAggoOD+dvf/ubtUJosSaaFEK6z2dD//gC9YweWSTfAuHHejkiIJun0in20LS2fCOpJ68qd0gDkBLeH0mQ6O+kk0LdB4hPnj40bN3o7hCZPxkwLIVyjNWnXTEXdfx+W9/8F48dT8MmX3o5KiCYpe6M5YDqrbRUDpksVRJgzehQelBk9hPAGSaaFEC7J+fxHWv/wsXPbjoVVs5Z5MSIhmi7HHnMqD1vXntWeZ4s0H0K0HZMZPerjfHqGTLhXfe+9JNNCCJekP/qys5xFGN05wOVJb3FS/n4L4XaBx8yeab+46pNp2pvJtDVVfhjrKiAggPT0dEmomyGtNenp6QQEBNS5DRkzLYSoVdGOfXQ9ZixfW4KV3uzmFO2hBP77X5g+3csBCtHEtM40e6ZbDqt+mEfhoAReWvAX0mhNUMshjG+I4JqgqKgojh075lx2u7koLCysVxLZVAQEBBAVFVXn+pJMCyFqdejV+ZTNKpoYeAUz32vPH/9obC9aJMm0EO5kL7IRVXzQuR01pke151oSLuRZLgRghKPa00QtfH19iY6O9nYYDS4xMZGBAwd6O4zzngzzEELUyufnH53l9Itv5IorzGOrVsGZM14ISogm6uSaZHwxljw+aelAiw4h1Z7burVZTkvzdGRCiKpIMi2EqJE9PYvo02sAcKDofu9ltG0Lw+PzSGAFD9j/wf7/JHo3SCGakJQV5hCPU6E1jJcGIiPNcjMboSBEo+G1ZFop1UkptUQptVsptVMp9ZC3YhFCVO/Qh79ixQ7ANusgBk0w/no/Gfg2KxjJP3gE9b/PvBmiEE1K7mbz4cOctjUn0y1bmuWMDLDbPRWVEKI63hwzbQMe1VpvUkqFAhuVUj9rrXd5MSYhxFkyF693lpN7XcKA0rfgASOHwDqj3GLvBi9EJkTTtN3ehyTupCf7KOxb80rGvr4w3/c2IkuO01qnkXl8La07BzVQpEII8GIyrbU+CZwsLecopXYDHQFJpoVoRF5p8QqbuZtBbGTS9b2d+ztdPQBKV6ftcGY3OBxgkZFjQtTXd4Xj+bF0Xo6vbq/9/NGOJbTFWLBl3/50SaaFaGCN4i+fUqorMBBY691IhBDlaQ0rVioOE82XTCLuJnOp4h7DW5OK8fRTkM4nf88Rb4UpRJOy1xwyTa/qZ8VzyvY3n0LMOSgDp4VoaF6fGk8pFQJ8CTystc6u4vg0YBpAZGQkiYmJDRugaDC5ublyfxuZI0eCSE0dCkCLFiWkpKyk/C0K8I8lsmgFAEvf/5rA6/tX2Y7c26ZN7q/7lJQoDh8eBSiU0hw/vpzU1JrnvAv2DXeWdy/fQHaPSn9K60zubdMm99c9vJpMK6V8MRLpT7TWX1V1jtb6A+ADgJiYGD1mzJiGC1A0qMTEROT+Ni6zZ5vlMWN8ufjiMRWO/9ghDg4ZyXSbNBuDqrl/cm+bNrm/7rNnjzFiCqBzZ8Wll46qtc66lu2hdHrK9tZgt94LubdNm9xf9/BaMq2UUsCHwG6t9d+9FYcQonqpi9YTjz+76MNFF1X+dWHv1hMOGWXb3oOVjgshzk3GV4l8zd/ZR08yIi4Grqy1ji3cHOZhPy2TTQvR0Lw5ZjoBuB24WCm1pfR1RW2VhBAN57JfHmcb/ckhlPHWxErH/Xt1NcsnDzVcYEI0UcWrN3IN3/Iof+dS2yKX6ujWMtm0EN7kzdk8VgDKW9cXQtSspMhBz9xNAARSSLdLulc6J7SfufxuWObhhgpNiCbLsr/c04c9ap5juoxqZU427ZOd6e6QhBC1aBSzeQghGp8Di/fTghwA0iyRRMRHVTqn7XAzmW5XcNiY/kMIUWehp8wFW4L6u5ZM+0aaDyBacyWZFqKheX02DyFE43R60UZiS8tHIgfRWlX+IKlj33CWMIZUIjlENA9mFREYEdCwgQrRhLTLNnum24yMcamOb9sIZ9k/T5JpIRqaJNNCiCrZ1m10lvNjB1V5jq8vTO2yhORkY/va0xAbUeWpQoha5J3Opb3jOAAlWOmY0NWleoHtzR+6wCJJpoVoaJJMCyGqFHHATKYDR1SdTANER+NMpg8dgtjYak8VQtTg6JL9zk+Djvp1p1uAa3+iA/rHMJU5ZBKBPbQj33ouRCFEFSSZFkJUYi9x0D17k3O708Sak+myOf8PyYQeQtRZxhpziEdaRC+6uVivRbfWzGUqAKH5HghMCFEjeQBRCFHJ4d8OEoaxilq6pTVtBnWq9twuXczysWOejkyIpqtou5lMF3RyYR3xUmFhUPZIQ04O2GzujkwIURNJpoUQlZz8zhzicbjVIPMvdRV6++7nRZ7hQ/7IoB9ebIjwhGiSrAfNZNrS2/Vk2mIxEuoyZ864MyohRG1kmIcQopKSNWYyndur+iEeAJ2sJ/kdLwOw69Bw4FlPhiZEkxWeaibTYYNdT6YBwsMhK8soZ2ZCq1bujEwIURNJpoUQlWzK6YlmLBewicCEC2o8N7xPB2c5Iu+4p0MToknSGh7Ub9OVnfRkH3eO63tO9T/MvI4ebCScLJI3rIYe51ZfCFF3kkwLISrQGl48dRdnuAuFg8N3O2o8v3U/M5lubTsJDofxubMQwmWnTkFi/lBgKGFh8FSfc6vfVp+mM0cBKDwp0+MJ0ZDkL54QooJDh8wxly1bWegUXfN77tadAsnAmOfWFxt5yWmeDlGIJmdvuVXEe/Wq8TGFKhUFmqsgFp7KclNUQghXSDIthKhgkzkjHhdcUPsfdaUg1dfsnU7fJkM9hDhXSUlmOca1hQ8rKA4xF24pSZGeaSEakiTTQogKzk6mXZEV3NEs7zrh5oiEaPr27yp2luuSTNtDzWTanibJtBANSZJpIUQFgz99hDd5iNv5fwyLce3j4vwws2e64ID0TAtxribOv4UUIllBAhfqVefeQLg5zINMSaabu99+g4svhn794JlnoLi49jqi7uQBRCGEk3Zoxh6ZSwTGH+PDXQ8C4TVXAmyR7aF0SXHb0VMejFCIpqld+k4iSSOSNPZ29zv3BiLMnml1RsZMN2cLFsCkScaz4ADbt8Pu3fDFF+c+Fl+4RnqmhRBOp1YfIkIbiXQmEXQe1dWleqptG3MjNdUDkQnRdBXnFNG5ZL9zO+qS3ufchjXSTKatOdIz3VydOAFTppiJdJmvvjKSaeEZkkwLIZyOLjQXazkQPgiLj2vdGD4dzGTampni9riEaMoOL07Cih2AI9ZogiKDz7kNvzbmJ0h+eZJMN1cvvADZ2Ub5xrbL+KTL0zzA24SRxYsvGlOfCveTYR5CCKfCVWYynd2z5pUPK+jfn1eYTgptCIyMZ5gHYhOiqUpduouy9Q5PtepL5zq0EdDe7JkOKJRkujk6eBA+/BBA8ypP8sTp1wG4FXiEvzN6+1JWrepCQoI3o2yaJJkWQjiFJpnJtN+FrifTQYP78jSvADDIATPdHpkQTVfRpp3OckG3uq1cGNjBTKaDimXMdHP0/vtgt8NdzOYJXq9wrCvJfMJkPp+/jIQEGZTgbvIVFUIAxsOHXTPMZLrjNa4n05GRZjlFRnkIcU4CD5rJtF//c1z6sKyNfj0ZyTL6s4VbQr5zV2jiPFFYCHPmQFtO8TqPmwdiY3FYfAAYwUpy5n0rQz08QJJpIQQAJ1dVfPiwy5hol+u2Kff8YUqKjMsT4ly0TTeT6dZj6tYzHdoumBWMZBv92Z1Xl4Ei4nz2+eeQng7PMYMwSgdN9+wJmzahH36EZT5juZaFzE27qsJqm8I9JJkWQgBwfME6Z3l/xGCXHz4ECA6GwECjXFQEubnujk6IpikntZAupTN5OFB0uezcZ/IA4+fPWjpws6jIeInm4+OPoQ2n+SNzzJ1vvw2Bgfi8Pos3rviNb7gWBz4sW+a9OJsqrybTSqk5SqkUpdQOb8YhhICileud5TO9hp5TXaXgLeuj/MrFbCeOzDVJtVcSQnB4cRI+GPOYHfONxi88qE7tKAUtWpjbZTM6iKYvLQ1+/RUe5G0CKH0XNXgwXHaZUbZYGDXKPF+Saffzds/0XGCCl2MQQgBhe81k2i9hyDnXH6zXczFLiGMnOUmypLgQrjix/jiF+AOQElm3IR5lWrQAH2yEk0l2hs0d4dVJQQFkZHjt8s3OggXGg4cbGMy6sEuMd1ZPPllhhZbRo83zly71QpBNnFeTaa31MkB+5ITwMq3hTjWHyXzMWzxIx+vPfXK7/BDzKcT8ZFm4RQhX/Gi5ghByiWM7myfOqFdbP5wagA1fMmlJ0c79tVdws+JiePrOFF4IfYP/tHqCB3v+wLq18gCFp82fb/y7kOtY/cJPcOQIXHddhXMGDjSG41mw0/boek7tl7F47iRT4wkhSE6GdRk9WEcPvmsxmQcuPPc2isPaQOlK4sXHZEoPIVyxeTPYsbKTOCIvqV9b2urrLBemNPw4j5cnbebhby+jDaVvpve/zrcXXcPSbz9l9BXnvhCNqF1KCixZYpSVMpYRp2NUpfOsVng/7EmuzPuACLLY8H/f0O6Fqxs22Cas0SfTSqlpwDSAyMhIEhMTvRuQ8Jjc3Fy5v17yyy9tAGNKrl69Mli2bNs5t5Hh62+W9yRVuJdyb5s2ub9143DA+vUjKPvrtRyFAAAgAElEQVRTXFKymsTEuj85GGIJdJZ3r91Kdmx+fUN0+d5u2BDB6W/XEkHFBWOudnzDomtvYt5HT9E+qqTe8YiKvvmmAw6HseRPfHwW+/ZtYd++qs8NDsomAmMO8swvvyJxbKj87LqL1tqrL6ArsMOVc3v16qVF07VkyRJvh9Bs3Xef1sZgD63/+te6tfHLxHecjayIu7vCMbm3TZvc37rZs8f8uWvTRmuHo37trYua6Gww8cEv3RKjq/d29Gjj0v3ZrHdFjtI5l1xn/udAz+nyvLbb3RKSKGfsWK07cExbsOl33qn53B8fWuS8H3sjhmit5We3NsAG7UJ+6u0HEIUQjcDJX3dhwQ7ARRfVrQ3fjuZk09YsGTMtRG32LdjBDXxBNAcZPEiXf16sTmxB5nQetoyGG+aRlGQ+1LbTOoDgDUsJ+ekrjt/+pPOc25Jf4qtnNjVYTM3BqVPG1/0bruEwXZmS9JQxtUc12l1jztLUOXOrMchduIW3p8abB6wGYpRSx5RSd3gzHiGao9xjWXy+J45MIviOKxk+1FGndgI6m8l0UI6MmRaiNv4L5vEFN3KQ7jyZ/Zd6t+cIMZNpR1bDJdOffWaWr7oKOpeuGdPxo5c5HJUAgC82Yl79I2mnvDfLSFPz5ZfQ27GDQWyiE8cImf0Pc7LxKsQmtOIA3QDwp5i81ec+nE9UzduzedyitW6vtfbVWkdprT/0ZjxCNEf7/7sSC5oW5NA18DQtwuv2ayGkSytnOagw3V3hCdFkhe81F0oKHh5f/wZDzWRan2m4ZPrbb83yLbeUO+DjQ/sf51KgjLHc8Y6t/HLz7AaLq6n73//gD/zX3HHttRAeXu35/v6QFGpOe5r24/pqzxXnRoZ5CNHM5S1e4Syf6jGyzu2Edmlplku8NOOl3U7Jv+eQNmIiqaOup+hfc6BEHnoSjY+tyE5s1hrndtTv6ji+qrxyq7aonIZJprMyNU9tnMQzvMgItZLLLqn4yZZ/3x4cvuVp5/ZFS19h1zbpna6vEydg1TIbt/GxufMPf6i13unOZjJtWy3JtLtIMi1EMxe2fbmz7DOm7sl0eDczmQ53ZBiPuTSk4mJSLpqI75/uoPXKr4lcvgD/e+/gdOwo9LHjDRuLELVI+nIHoRhz/Z6ydKDtkM71btMnwkymLXkNk0xvmpfEDXzJizzH95arCGtR+ec+9j+PcSKgG/O4mQRW8uiTjX4isUbvf/+D8fxM+7L5SNu1g0svrbWerd8FznLA3u2eCq/ZkWRaiGasKKuAnllm70S32xPq3FZIZCB3WuYwkQWM41cK8hs2mT429VnarPuu0v62B9eQEjcWfep0g8YjRE1OfbXSWU7ueBH1fvoQsLY0k2lrfsMk0xlfLnGWj0SPBh+fSueowADSf97EZDWPY3Tixx9h0aIGCa/JmjfvrCEekyfXOF66TMjwOGe5VcouY35GUW+STAvRjO3+YDn+GE90H/CNpfOQtnVuSyn4LnIqXzORFYwkI6vhfr0U79hL+0/fcG7PCbiHN1s8hw3jD3vbM/s40e8yyMtrsJiEqIl13Spn2Tas7m9iy/NtZSbTfgUNk0yHbzGTacvFY6s9L35EGHfeaW4//DAU1X1K7WbtwAHYuy6TiSw0d7owxAMgemgkKRir1QbY8wk4LZ0M7iDJtBDNWM5XPznLR2Jr/4iwNi3NkR5kNOCw6YN/eg0fjB6WFT6jGbf7Pe4+9QJvjfjCmVC3Sd3JD39ZUVMzQjQIraHrcbNnut31bhgvDdjHjKMjxwglm/s7f+OWNmu8nh26ZW50bre7seZhYi+9BGFhRnnfPnjntQJPhtdkffYZ3MR8Aih9NzJoEMS79gBrTAzspC85hLBeDUXlSAeDO0gyLUQz1m6bmUz7XX1ZvdtrZU7oQXoDTeihNdyV8jIP8DaruJCjd79El66KwEB4aMlE/nPBvzhOB8aQyMR/XcbatQ0TlxDVOfjzAbo4DgOQRxDR1w1wS7uhbYM4QUdyCeVMjuf/vB/YdIZu+iAAJVhpObJvjee3aQMzZkAIObzG49z01xhO7Gn4Zc/Pd/PmwRTmmjumTnW5bkQE3B25gBZkM1Sv5WBoXO2VRK0kmRaimcrYcYKeBcYDKMX40uee0fVus0LPdHrDjJnesgVW7GvLuzzApcGruGrWCOcxqxV+v/wuborbxSoSKC6GG26AFJkGW3jRkf+Yb2J3tR2LJcDPLe2Wm8yD7AbIUY9+u8Ush/Y15l6rxb33aFYFjuNx3qCTPsrWa5/zZIhNzpYtYN+5m+EYvQLazw9uvvmc2mjfOxwwxugfORLk7hCbJUmmhWimdv39R2d5e4sRREQF17vNKUdf5BBdySaUNgv+Xe/2XPHpp2Z54kQIDa14PCgI/t/XYUREGNvHjxt/e2wyO5fwEv+li53lolH1H15VpnwyfeaM25qtVv7Kzc5yZvRAl+pYfRU+jz3i3L5s79usmrXM7bE1VbNnQwpteIJXOdaiN+qaayp+JOiCXr3M8okTgW6OsHmSZFqIZuqXPVEs4nKK8SVn5BVuaTPMmkdXkgklF3tawwyaXrDALE+eXPU53boZH42WTZiwf8kRlo98Sp5kFw0uPx9mZt7LO9xPEr2I/lP9h1eVCQiADj6n6cE+4oo3UlTo2U+HAveYybRloOtDVfq8cBNb2xv/bwuajs9MIe1wrtvja2ry8+GTTyCDVrzOE+z9aifMmXPO7URHm+WTJwPcGGHzJcm0EM1QXh68uvlSrmQRbTlNpxfurL2SK1qZ4zxUAwyaPrb+JPrAAUATFATjxlV/7mWXwfPPw+UsYjMDGbtmFjv+8LrHYxSivCVL4PuSS3mQd7ihbxIdL45xW9tKwUF7Z/bRi40MJjul0G1tV6VN2i5nOXx0f9crKkXnnz4kSxmr9XWxH2LrmAfRjgaem/488/nn5icO3bvDmLGq8kdxLujWDfqzhVv5hMvWvgMyo0e9STItRDP0ww9QWPp3tkOfCLoPqn4J2nNhiTQ/bvQ54/me6dMzP+QAPThENDO6fYRfLUNPn3kG7oheQiuM2GI//gtH/y/R43EKUeabcpNsXOGeD4QqyLWYYz3yTnpu4HRRoaZL8V7ndvuxsedUPyKuI4cfece5PS75I3753Qdui68pev99s3znnWCpYwYXHQ3vcR+fcBt3HHkdtjf84i3HjsFf/wpXXmkMz5s5E5KTGzwMt5FkWohmaO5cs3zDDe5r16+t2TPtm+P5ZNp/tTHHbVeSiYmr/SEuiwXGrp7JBn9jKjIrdkKnXE/a8t0ejVMIgOJi+OILc/uaa9x/jXxrwyTTBw9obmI+D/IWc0IfIqDLuc9RP+D1yazraY7NGv3lA6ydtaSGGs3XypWg16xhFEvx9YUpU+reVrdusI+ezm29d1/9AzwH//sf9OljzOyyaBF8/TX85S/Qsyfcf3/DPDzrbpJMC9HMnNybw6+LzNUSXJzr3yUBHcxkOiDfw8M8HA66pJirN3ae7NpS6C3b+hKwcD6naQNAuCOT4vGXk5N0wiNhClFm2X8P4ZdxEoDOneEi90wvXUGBr5lMF6V6LivZu9/CYibwDg/y2fA367aCo1IMWPcBe4ON8dZ+lNDnqWvY/O91bo72/Pfaq5p/8GeWMobtrUbTLmNX7ZWq0bo1JPuayXThjv3uCNEl338Pt94KOTngRxFX8w0v8Bwv8gw3lnzCvPfSiYuDX39tsJDcQpJpIZqZA396naO6I6/zGDcPO0T37u5rO7izOcwjuMizPdMZa/cRqnMAOE0bel/ayeW6cROi2PPG9+RizGDSoTiZzAFjyNp2xCOxCgHgN+MvHKUTX3EdD4/fXueP6WtS5Gcm08Vpnkum95XrzCw/O8S58gsPIiJxIad8OgAQSi5H7nuVX36pZ4BNyJYtUPDtz1zIGgB6pq+pOHXLOVIKctqZyXTRjobpmU5Lg9//3ljsZyy/sd8ayzdcW5pKv8wn3MYJOnDB0YVccgk88YTxac75QJJpIZqRouwiui/9kNak8xh/4+4hG2uvdA5Cu5g906Elnk2mjy40Y98fPhhfv3PrGRv96GBWPPi5c4XEzoX7KBg8gpO/1r3HR4jqHNucyrBjX2LFznUs5NqrPDOTTHGAmWSVZHiwZ9ocLl2vZBogcnAXir79mQzViq304/f2OVx+uTFGWMsziTz1uI2/8ahz23LnHRAVVa82bd3MZNpyoGGS6enTjZVxr+YbfmQCnWyHK51zRoXzC+PRGl5/3fj0pvz3WmMlybQQzci6+/5Le4cxnOG0pR3DZ17r1vbDu5nJdIQj3aNP5xes2OAsZ/cYVKc2Jrx1OT/d9QVFGOOt25ccpcX4oWyZ/plbYhSizJ5738Yfo5ttT4shdLvuHGa/OAclgWYybfdgMr1/n/mz3aNH/dvrcnkfsr76jTvbLyKbMGw2uOceuP56SE2tf/vnqx9+gG6//Jt4dgDgCAw2BhjXU0Bf86YFnT5odBd7UHKy8axOD/bxKbfiR4lxoGVLePRR4/80bBgBLz/HyAkhznobN8KYAVl89EFJo35jJcm0EM1E0ZlCoj+b6dzeffmjBIT6uvUaQS0DyMNYUcsXG/kpnps7NiRpk7PsnzC4zu1c8cFEVk7/zhl3MHm8+2ou99xj9KIIUV8pe7MYvMacuSL3rkdqOLt+bMFmMu3I8lwy/bfVF7KPHvzCOHr6HHRLm90m9mPh+o4MLLf+y8KFmo0drmbJdW9TlF1UfeUm6MwZmHnHAWYx3bnP8szT0LFjvdvuEBPKKYyHRq32Yjji2SFub70FdrtmLlMIIc/Y2a2bMYbljTfgpZdgzRpaTL+X77+Hf/wD5+xMbxTcy7C7+/PSqJ/IzPRomHUmybQQzcTaG98gymbMPZSuWjFo9p88cp0zFrN3+swhD2WjWtMx05zOqf2E+vXyXfzKJez5aA0HfHrxEVP4kDt4/32IiTE+aszJqW/AojnbNmkG4RgTBCf79+SCV2702LUcIWYyrc94Jpl2OKBz0T56cIBx/EbHHu5bRa9jR1ixAu67z9i+mc+YYPuOsQsfIjMimt8mvEr6Xs/PYe9tWsPD9xTx2snbCMXolLD1iIFH3PNGLDq64oweFQbBu1lREXz0EUxkIQmsMnb6+hpT23Q661kXpbBY4OGHYd06mNxlBbcyjz7s5tkVl7G+/TX876W9jW4FW0mmhWgGDny9g6E/v+Tc3jbpRULbh9RQo+4ei/6SfmwliqOk+NVvXF918g6nEuEwEvUcQug2pnO92xw0JZ4WSev5/vL3AGP8dVqa8RDMc+0+YGXcNJLmrpaFJcQ52fyfjYzebvZKpz84A4uvj+cuGGom0yrXM8n06QO5zrnai/AjKPrcp8WrSVAQvPuuMcThoUBz7ul2jpNcvHg6LWLasa7t1Sy/51My9jXNxPq1WQ4unTfF+dChw8eKdd7HxjKXbtCpExyg3NPnhw65pd2q/PADZGXBNMrNI/7gg1T4CKIK/fvDh88lU+hr/q26tOhbJj0bS2LEdXz/5DIK8hvH72NJpoVo4s4kZ8JNvyMA4yPSXYEXMPL/pnnseieihrKdfhwniowznkkajv1kPiR4KKA3/gF1mJarCpHdW/D590EsWABduhj7/CnksfwXSNg5m5ipF3HKrzMr4+5mwzMLyUxKcct1RdN0encG4X+6GV+MbrQdLUcxcNZNHr2mahFKIf6k0YoCm3uHcZU5vd4cEpDi36nuq4fUYsIEuOD0jyy75nVOW9o59/tiY2jKd4x8fzLhvSLZFTSYJUOfZPHbSRw+fH4/tKg1vP5yMVFP384tmM9uqFmzYHDdh7OdrXNnOIS5rrg+6Llk+pNPjH+v5ls+nrQQrr7aeBrRBf5/nEzAkX0kXzzVuc+CZnzuQq58bTQnQmNYNOhZ1nywDVuJ9268JNNCNGEZe9M4Gn8F3YuMRUkKCMD3s4+x+nuuZ6yVOTsenlpR/OjeAnYTiw0f0tv0cWvbShkrciUlwezZcF+bL+iIOQd1e/sxEnZ+wOCXryMiti1H/bqxOvoWlk2YybqnFrAvydHoPoIUDe/UjjRODbqCaLsxh28uIbRa+CHK4p43ftU5dNUDBFJIJGnMi53hkWuUn0IyK7T+nwrVxC/Un1FfP0ZE1mFW3TmHHcHDKhy3oOlTsJGx61/jlYdOEh1t/A4aOxbunqb5beJbrH32Ow4sSiI/Ld+jsdZXaircfDM89YyFUMyxZfZp96Aede84+4gI2Ol3Ad9zBe9yH3kDEtzafpnCQmNhFgA7VgY+f62xDGjr1q430q4dXX6dQ8Gy9RzocVmFQ90d+7hi00sMv7s/D4XN5corYdYsWLrU+GSxoVgb7lKVKaUmAG8BPsB/tNazvBmPEE3J8uWQeP2nPJuzxrlvwz1zGHlNb49et6U5ZNpjyfTP1st5jcvxpZgXfpfLWA9cw9/fWLLXMXUyG97rRu7bc+h/4EsiyKpwXqeSQ3Q6fAgOf0b64pa0njURq9UY+9mpEySEbOXy9I9RnaKwto8koGMrgjq3JjS6NRE9WxPUOqhuC16IRklr+O472HLrHJ4tWAuAA8X+p+cwYKQbpr2oRWgL83sp10PP/xbuNdd9zm/TxTMXOYtfqD8XzZ4Ks6dyZMkBDr8yj4jVi+iduw4rxkwUmzGGDWRmQmIi7ExM5d88DF+b7WSpcFL9OnImpCMFER2xR7aDVq3waduazGv+QMuWmK8IjZ+/Z382HQ7YsraITz7z4YM51tJ7ZuVmPmN9i3H0nNQfv3++4/bfEUrBzm5Xc9WeqwEY1Rf6ufUKhuXLIb/0PUzPntC3b93bChw5mO77fuTM6l0kP/Im3dbOI0Sb3+SLC0ZyYJGZvFspYbdPHNmhUeS1iaYkKhpL966E9GhPREwbwnpEEt69FdaA+qfCXkumlVI+wHvAJcAxYL1S6huttUzyKkQdZWU4WPyzhblz4ccfwco9XM0c+rGNpTe8w9h/3uLxGFq2BB9shHGG/BM+QLjbr7Gr9LdECX50G9yy5pPryeKjGPzgRfDgRRTlvc/GD1Zx5tPvabVnBTG5G53DZwD2EAsobDZjKqjkZIhlHaN5A9ZX3X4BAWRbwsn3CWVVxJXM7v0PQkIgNBRCQmBw+mJi0lehAgMgIAAVGIAl0B9LUACWoAB8ggPwCfLHGhIAnTph6doZf3/jzYDVCr75Z7BiwxroizXAarwCfT3eQ9qcaIcmaVMeS9aH8N//wtq1oHiMeFZxDd+wYep7DH3Zcw8dlhdS7lEIjz04m2z2TOsoz/ZMV6Xz2O50HvsM8AzZR8+Q9P4S0lbsYYhvGBs3GuNzAXpS+aG6cJ1FeFEWFO2EdKB08b90WtL6w4rLwV7OD3zF9eSpEPItIRRaQyiyhlDkF0KJfwg2/xDsgSHogECyWndn4/D78fPD+WqXvpOOR1aDnx/2EgeOwmIcRcU48gpxnE7BN+UEYWn76WffzD4+Jxdzffnb7gqi28xf8GvluTfbnTvDnj1G+cgR6OeBbLossQW4/HL3tBl2YR/6rf4Ae+5b7PzHDxT9v/n4Je/lQEnFFchiSKKHfS9k7YUsYC/wW8W2HCjSVCuO+vXgngGradXK6DRv2bLip6y18WbP9FBgv9b6IIBS6jPgWqDaZNqWks+qB8+a/7WawVFHhk5CW32dp1jzs4na/G2l86qqrn2sJA+/qcLxgKxTtN/xc8XzqomzJKAFRy+oOH9vSMpB2uxdXnPFUgUt2nI8fkKFfeHHdtD60PpK8apyjZUdy4nsxqnYMRXOa31gLRHHd1S+WBVfgIwOcaR2H17hlA57fiM09UDFa2td5dcvpdtwMqL6VWi+y9ZvCDpz1nLNZ9VNPX2aZR/s4XjsOLIie1Y41nP9p/jlZ1VVrZJD8deSF15x6qC+K/6NxV5xnkp1Vktlx5KG3EZRUIRzv8VWTNyK9ytfqJrvvW0J9+Cw+jm3/fIy6bP+v1XWP7sFh8WXrSPvr9B8SOZRem37vHL9omJIT8fnTAYBWadol7mLZbaLuI1PnKfY8OWhwNm8cu9Rxr5xfZXxutul2//GqzwGQOLix2HGa26/RvlJ/GNj3d58tfyDrQz68yj48ygACrOL2fb5VjJ+3ojes4fd+V3pmA/Hj5t1OnG0xjYDKSTQcQocp1iVksLSs4Zhx/ITY/i7S/G9zmM8wesV9s3nLn5H5e8fOxZsWCnBFztWbMrKKyEz+azFNKxWnK83Tk4mpmgbDmVBKwsa41/ntrLQwu5gi68fc6NnsL3laCwWnK8/7/gj4cWpznO1xQJVlLFY+Kb/s6SF90ApM3+4feXdWLQdUKXPhSp02UGlKHtYFKX4dthLFAS2dB7yK87lmjVPm+eW1i9rXJfVL93//ahXwcfHeXpo7klGbHrLWV/ZbKiCfCxF+fgUFeBTmEuLnOO0KzzMbsZwLwucX1+NhT+3+j+6PL6aoU9e6tL9c4fQULPsqWTa/7SZTPt2b/hkurwWncIY8vJEAC7H+L156JDxhjtjeTjLvvsToSf3Epl9gDb2E+b8xmfJoPKb8hByCaCIAF1EK3s62IEiKJvZrbwVJPDSL/dX2PcAv/I2D7n0/5jMJ3zLNcTFwauvwhVXAKUrtHpK+Yk0PDUz3uLF8CBvcYhorhp/KeCeBygBfEIC6fvs9fDs9WgNe/cbPeErVsDWrdBv+1aqud1OFjStdRo5RSGsXVuPYLTWXnkBkzCGdpRt3w68W1OdQcbPiUuvELIr7OrFHpfrZhNSafdolrhcfw+9Ku2+hU9crr+E0ZV2P8IbLtf/hFsq7X6DR1yu/waPVNr9KTe7XP8R3qi0ewmjXa5/M59W2p1ET5frj2ZJpd3ZhLhcvydJFXaFkO1yXY13v/e2Eu/cVErr3/1O66NHdYNKvPV9ZzzLYu7UWmu9ZMkSt7Vvs2nt62v+t3Ny3Na02+TlaZ2UpPUvv2j9/dMr9G/jXtKJfe/VKzrdpDdEjNO7A/rr45YoXYB/hfv3L+6udFv/zV0u3/8X+Uul3V8x0eX69/FOpd2rGeZy/ev4stLuZDq5XH8YqyvtLsHH5fqdSK6wK5LT5/Sza6W4wq5+bHG57mkiNTg0GN+f996r9enTDf+9t3NDvr6QlfoSFus/dv653u1V9bO7MWSk8/+9+536X6Oh2EvsOmXHab1n3ia97rlv9bLJ7+vfxs7QiQMf1l/2e15feaXWF16odUyM1pGRWt+hPnT5/v/K2Eq7H+M1l+snh8frH7+3aZut4b4eL7xghjB9uvvbT0nROpgcnU+A1qAdISFanzrl/gtVw5ZXqI/8sEOv++t3eulN7+plQx/VaztP0lvCRur9vrE6XbV0fgE209/5tfClSN/P2/prrtbABq1rz2m92TNd1ecWutJJSk0DpgHUbY0zIZqP7hwgtmcWwxMyGT8+hY4dC9i/H/bvb7gYUgrLPTiTfprExERyc3NJTEx0S/tZe3K5t2Q/B+nGyRbd2LChAZ8yOUc+PhB0CXBJAhqjk6QEyAFOAdqhKTlTQlFKAfasYkJ0MH/330JBgQ/5+cbLljSEhceDUMUl+JQU4VNSjE9JMVZbkfGyl+BrL8LXXkRBaDu6BOZRUmKhpEThcCiKswPJKInAig1fSkr/rfoJyRIqz/5gwfVlrx1VPNN+LvV1lX8WXOfN+v4Ucc3gnXQe7s+4cSmEh5ewa5c5JKmh5G5NZRW/AyD5WFcSEz+qX3tV/Ox2yzfHTB+0pXPKTT/bDaYd0C4ExsagiEEDLYHHSKxwmnZE83PeT5RkFmM/U4wtqwhHdhGO7EJ0TiHkFWLJLYDCEjID2zK18yFsNkVJiQWbTRF6rB2Lj9+C1VYMPhbsVl+01Yr29cUREYqjfUuCe4XhO7QLtvAw/FnO8uUN92XIy2vHDexgLEsYPHsPW9peTdaAAW5rf/ny1lzMQQIpNK4XGcmG3bth9263XaNWAcCYYBjTF+hLPpAPZAJHAUexneITeeSn23nbbxPZ2b5kZ/kwcE0aK1vcD4sqj2ioijeT6WNA+dm6o4ATZ5+ktf4AjMkJY3zD9ar2E84+xfmRX/lfg9cNtVJS7n8XUdiClZtvdSmwEp8AbjVHOaAUdDjTlhU7bqt8chW/e88EtOO2QRWHOfVIjWbF3t+7dP3UsFj+cNbYpQ7H+7L80JRqrqkqFEtaD2PqWR99Bx0cyvITf6zyemd//YI7DuOObhXPKUkay7LUquYlVpXi6RQdz13lphdWClK3Xc3SM5U/jy9fNS8/j+CgYPrHdiesTcWv3971t3KywPz8u6YhZGP6daBvRMV961ZOw2ovqnzRKm7gpCFh5ASZ21abL4mrHqj63V+5QMpK0xJ8sZX73gvJj2Dphqo/6tOoCv8Xm48/D4+oeE7kmSgSt/65cmUfK7RsiSWyFb7tWtFqeE+6XBrD7lA/jHHK0ZXrNIANq4rhK6PcwpHPBWPGkJiYyJgxY9zS/pbVP/AmxtdjI2MZNOa3Wmqc72JcPnMA8HKlvZ9W2qMdGnuJA1tBCbZCG/YiG7aCEl7wC+ZZP7DZoKTEeFmT57E3Px+HzQEOB9puvBw2s7x3zx56dOvBY1368kCo8VBV2evYprmcKCyoUFc7HFBadu53OHg4rhsFIUY9MPqJVi7/N0o70NoYmqV1aZ+LQ5dul54IPDMknJIA5yY+RSEsXf0WaF25vtbOoWpl7bwy0gIWs35gdjsSN77i3KF8fFDBQajgICwhQfi2CCS8d3siB3UmvFcbvvYpezNRcZhaQ8psd9pZDta59f65O/tn11ZkB4c5jmn81GvwC3Pfoi3ns7sq7bm99NU42e1w4LW/M43Zxvhxn4ngpt/TAN9+C+PLzS0dcuONbvs74HkvcDPwmovvr72ZTK8HeiqlooHjwM1Azdlut7ZclPRZjaeUuajSnvZQbixpbcZU2iauiqUAACAASURBVNMb+D+X619Zac+FpS/XVH5UZULpq3YjgT9U2ntT6cu1+pU5PyCo1agq9z5aa73ExERGjRlTTf3nXbp29df/m8v1R1faEwC8XY/6bYA3Xa4/rtKeHuDimNnGILC9+cBhQGFWDWfWTd52c+ninMhuNZwpqqMsCqu/j2tTJPbpXuspGYklDBhT+TsfgMvHn2N0Z7nrDpdPrRxBEPCgy/XHVNrTFnBtPtzGIqS9OWg6WOcabxbc+PxaSpqFISTTnpP0aJHKZ5JIn7c6d4Zfyne6uHnhluXL4Y/8au4YX8/fBY2Y15JprbVNKXU/sBhjarw5Wuud3opHCOEewVHmxwLBxZlub9+xz3wQ1tal9kRPiObEt0Ugdiz44CCQQgpybQSGuu9P/clTihN05AQdKenqtmaFF0RFVV64xV3vu/Ly4OTGE/QtnVNC+/mhEjwzl3Vj4NV5prXWi4BFtZ4ohDhvhESZPdOhdvf3TPufMHum/WKlZ1qICpQiT4XQQhtLieeeyiUw1H3TU54oNxizQwe3NSu8IDAQ0kM6QulUzSX7DuFXcxWXbdgAox3mEDyVkGCsE99EyQqIQgi3Cuti/uEO01nGWFs3aplh9kyHXSA900KcLc9iDvXIO+Xe+fFOnjTL7du7tWnhBbmR5jsiy5HDbmt30yYYzVJzx7jKAxibklqTaaVUkFLqWaXU7NLtnkqpqzwfmhDifOQb5Etu6fyoPjjIOem+Zdi0Q9O+0OyZbp8gPdNCnK3Qaj4sXpjm3mUQs/an0ZJ0QEvPdBNg7RBKHkaPsTX3jLF8pBts3gwXstrcMWJE9Sc3Aa70TH+EMU152dNzx4CXPBaREOK8l+Nj9k7nHHHfuOnMpBSCMdamzSKMyBjPrn4oxPmo0Gr2TBemurdnesiiF0inNUX4M+7gbLe2LRpeZJti/n979x1nd1Xnf/x1prc7M5maCZlkUicJhBQiRSkTQBRUEBGxoCL+jLoqrmUtiyvsur+ffdWVtWBF18aiAq6K0gaQGgwJCaSQMunJ9HKnz9zz++PezPeGTLm9fO/7+XjMI+d77/ee+cDJnfnk3M85p4UG54EYLULctbGHFSfqpbOzYd26mPSbqkJJphdZa79C4BwZa+0gk+8RLSICgDfHvwhxjGy8R3pj1m/rJuc0weP58+N1yq5IWhvJc2amRzpim0zndfrrPPIYpTho5xBJT7W1wzFPpgcGoHzn02QFjg7xrVwFxfE9zTHZQkmmR4wxhQQOVDHGLMI/Uy0iMqlPrnkQD73kMsrRqpUx67f3hUMT7e7S+mnuFMlcIwVOkjvaFdsyj5Jep2i6eJGKptNddfXQSTt60NISdZ9bt8I2u4IPcRv3lr6D7GuujrrPVBfKbh63APcB9caYXwCvAm6IZ1AikuZqak4sEKe7G2bNmvbukA3vdmamhyrnTnOnSObqKZvPjqON9OGhbzy2OyiUDR2baM9aPjumfUvi1dYOszfGe01v2gSHmct3+BBdr/sQV34u6i5T3ozJtLX2fmPMJuBc/OUdH7XWpu75vSKSdMHJc1dX7JLp3Symheup5yD9i1bFplMRl7nn0m9z2w5/+5vzQj3ua2bjY5aacWdmumqlZqbTXU3NMD+miZv5dzo9DXz3H9ZE3edzzznttWuj7i4tTJlMG2Ne/r/gxDtonjFmnrV2U/zCEpF0Vh60rW13DLeafjDvcn7B5QD86E2x61fETTxBpczeGFZ5dB3oo4pBAPoporhSNdPprqJimC3Z57Fp/Czog28s9J/5G41NQdlhxifTOOcvFwDrgC34Z6bPBJ4G3L3PiYhErKbIy2KOUk437KmAGP1APeSUTFOvkmmRSZU46w/pi+H6w+7tR6kKtDtyZ1OsFcBpLzvbf/jOwUAF3aFDsHhx5P35fPDiC5YT+1SceWb0MaaDKRcgWmvXW2vXA/uBtdbaddbas4A1wO5EBSgi6ef8HT/kJZaykbNZ+9i3YtbvQadkWsm0yBSCZ6ZjmUx7X3JKPLoLVOLhFvPmOe3gn7GROHAA/nHoi+xlAX/Mu5qqjX+OrsM0EcoCxGXW2q0nLqy124wxq+MYk4ikuexKp84jp6+LsRj06fOdPDM9V+sPRSZVN9zC23kcD3007FoIXBaTfof2O4sPvR4tPnSL4ImJAwcAa4l039Ht22E1m1lACwtGWuBIZpzxF0oyvd0Y80Pgv/Fvj3c9sD2uUYlIWsutcVYc5g50xySZ7tjVwY9HbuIQczlY2EhJyY0x6FXEfeYdfJxfcD0Aj+14K7FKpscOt060R8pqYtKnJF99PXyKL/NWfs3yD+yDgu/DdddF1NeLL8Lr2eI8sCozFoqHkky/B/gg8NHA9aPAd+MWkYikvYLZzsx0wWBXYMlSdNqfbeEd/BKAXXYloGRaZDI5s5w6j5yh2K1AHO4aYIxschjHV1kds34luerrYYxjrGEzDBHV9nh7nu9nCS8B4DNZZJ1+eoyiTG0zHtpirR2y1n7DWnt14Osb1tqhRAQnIumpcI4zM100EpvtPHqCDmzp8ajGQ2QqeZXOCsT84dgVTd+97DPkMUIFHey84mMx61eSa948Tj4FMYqDW0af2zZx8uFAfSMUFkYXXJqYcWbaGLOPwOmHway1C+MSkYikPU+9MzNdMtoVkz6DD2wZrNTqQ5GpFFQ5M9P5o7GbmW5rA0sWXVRQOm/m+yU91Ndz8imIEc5MWwslezZPXJvVmVHiAaGVeawLahcA1wIV8QlHRNygdL4zM13qi83MtO+gMzPtm6OZaZGpFFQ7yXTBWOxmptvanHa1qjxco74+NjPTra2wZMCply46N3OS6VDKPDqCvg5ba78JXJyA2EQkTRVWlzBGNgDFDDDaPx51n7nHnGQ6p0HJtMhUimqcMo/C8djNTLc66w+VTLtIVRUcy2+YuLYtLf7tk8K0fTusClp8aNZkzsZvoZR5BB+3kIV/plrHHonIlEyWoceUU2k7ABg+Hv0yi+JuJ5kublQyLTKV4tnOr+gSG7uZ6aWHHyKXctqporpyLiHMx0kaMAbK6kvp2F1BJZ2YkRE4dsx/mksYtr/g43qedx7IkJ08ILQyj68HtceAfcBb4hOOiLiFN6ecylF/Mj3SGv1+HmUDzh63pY06MEJkKnmziifaJXgZHrLkF0R3WqFvzMdvui4jB/+nTMOzhoG8qPqU1DFvHuzbvYBKOv0P7NsXdjLd9sw+PPg/CekvrqZ4dubsRR7KPyvfe+I0RGvtq621G4CReAcmIumtJ7+GNqp4icUM9Yb/keHLVYw6yXTFisz5IS0SLpObwwD+XRSysPQdH4i6z5793ROJdA+l5HuUSLtJLOqmh5/fOdHuX7wq4oNf0lEoyfRdIT4mIjLh0+c/QQ1tLOUl9nuWRdXXQOcQ5fQAMEY25Qu1BlpkOgNZTt30YGv0pR5du5zVh905Kph2m1js6HFH2xWU0c0FPMrQp2+NXXBpYMoyD2PMMuB0oMwY86agp0rx7+ohIjKlWc6GHni9uVH11f5iKyd24mrPqmV2tmo1RabzTMnF0NuLlxLO6I/+/eLd5yTTvQVKpt2mvh62RDEz3dsLhw8DlPFM3gXMuTaGwaWB6WqmG4HXA+XAG4Ie7wPeF803NcZcC9wKLAfOttY+G01/IpJ6yp2tpvF6Q1meMbXj/SV8jW8xm2NUzs7n/VHGJuJ2/7b81zz9tL/9eAwqMgYOtE+0B4urou9QUsq8efAFrmQ7y6k7bwG/+G54i7y3b3faS5ZATnQ/8tPOlP+51tp7gHuMMedZa5+M8ffdBrwJ+H6M+xWRFBE8M93XF91P1sODFXybmwB4/VqUTIvMwBO055Y3BrvjjRx2ZqZHSjUz7Tb19XCIeg5Rz+I2IMwPE4OT6eXLYxpaWpiuzONT1tqvAG83xrzt5c9ba2+K9Jtaa7cHvkekXYhIipubdYRXs41ZdFHbUgjB9XhhOuasPaS2NvrYRNyuxCmZpi8Gu+ONH3OS6bEKJdNuUx90qOyhQ/7TDMNJ0Y48fZAmdrOFVSxfnnlrWqabLjrx7wyVYIhI2Fbuu5cP8kEA/vzSO4ArI+4rOJnOoN2WRCIWPDMdi2TadDhlHqZKZR5uU1rq/+rthaEhaG8P72Ceqsd+z8N8FICdz9wEfCs+gaao6co8/hD4845IOjbGPABM9mvv5kAJSaj9bAA2AFRXV9Pc3BxJOJIGvF6vxtdF2oec3+B5A71Rje1zzy0BTgOgr+8lmpsPRxmdxJLeu6lnwe6j3MwePPTR9ae1NDdE9pHOibEdPXpw4rFWO6rxdong925l5Tp6e0soZICHv/NH5qzOY6ysLKR+Svf9faLdVZadcX8/pivz+ANgp3reWjvtNJO19tIo4gru53bgdoDGxkbb1NQUi24lBTU3N6PxdY/nnhqG3/rbxSN9nBvF2Pb+n4/xcTZxnFpmzfm0/p6kGL13U0/u8E28im8DcP/AN2lqui6ifk6M7bOjX5p4rGHdmbxC4+0Kwe/dZcvgM/s2sIEf+LeI+MlP4KqrZuxjaAi2DXxy4nrNDa8nP8P+fkxX5vG1hEUhIq5TWOds51E02htVX/Wtf2cNjwHwfP4HoupLJCPEuGj6qKljN4uopo2CetVMu1F9PXRQ6TwQ4vZ4L20f4wy2TVznn505x4ifMF2ZxyMn2saYPGAZ/pnqndbaqE5ANMZcDXwbqAb+aIzZbK19TTR9ikhqKZ7rbOfhGeuOqq/yQadouqxRRdMiM4rxdh4fL/8JuwNl0zsuiLo7SUGnHNyyZ09Irzv88C5WMgxAW8FcqisrZ3iF+8y4X5Ux5nXA94A9gAEWGGPeb639c6Tf1Fr7e+D3kb5eRFKfp96ZmS4d74m4H2uhYuz4xHXFcm3nITKTrFJnZjprIPqZ6c5Opx287aW4x7x58ABLnQd27Qrpdf1PbJlot81ZRSZ+bhHK5q9fB9Zba3cDGGMWAX8EIk6mRcT9Suc5yXQZ3YyPWbJzwt8Os691kDL8ZSKj5FBSr9/kIjPJLndmprMHo5uZ9vmgO+jDJSXT7lRfD7tenkyHsEdezgtOMj26PPNKPABCOWO09UQiHbAXaI1TPCLiElkFefRTBEAO43QfiuwXevsLzqx0e3YtRkeJi8woZ5aTTOcORjcz3dvrT6jBXz2SG+aBHpIe6uvhKHX0EfhUo7sb2tqmfxFQfei5iXbBeavjFV5KC2Vm+gVjzJ+AO/HXTF8LbDTGvAnAWvu7OMYnImmsL7uc4vEBALwHu6hs8MzwilP1vuQk0z0FtdTFLDoR98qtcMo8coejm5nu3tXK+7ibTioYLaoHzokyOklFc+f6D9PbZZdyFpv8D+7aBTU1U75mfMyy2Osk03WXr4l3mCkplCmeAuA4cBHQBLQBFcAbgNfHLTIRSXveXOfzYO+hyBYhevc6yXR/ieqlRUKRX+n8wzV/NLqZ6eEtO7id93MX13Kr9xPRhiYpqqAATjvtZaUeO3dO+5oDTx2hBv/sdZ/xULp6YTxDTFkzzkxba9+TiEBExH32l51Jx1Ax3ZTj6Q/lg7BTjex3dvIYLtdOHiKhKKhyZqYLokymBw87qw8HCzLvqOhMsnAh7DzU6DwwwyLEPS8MsZFrWcNzDJTNYVVWZpbhhbKbxwLgI0BD8P0zHdoiIvKdV/2S3wUKwe4M7SCtU/iOOjPTvmrNTIuEorDGmZkuHI+uzGP0uJNMD5comXazhQth16Ohz0xv7FzEP3MnAB99xxjfjGdwKSyUqaK7gR8BfwB88Q1HRNwkeNV/V1dkfWS1OTPTZraSaZFQFNdX8DPeSR8e2rNmc0sUfY21Osn0uEdbebjZokXwe/wz0z0lcyibYeuWF1902svOiOzTRzcI5b98yFr7n3GPRERcJxbJ9K9rPso3t7+GWo7z1vPPi01gIi5XOLuM92T9zL8LxxjcPAY5EeY6vg4nmfaVa2bazRYuhC2swkMvl13m4bc/mf7+F15w2qefHt/YUlkob61vGWNuAf4KgSNuAGvtprhFJSKuUO5sNX3SPrXheK5/Kc8GFsS8V5sIiITEGP+J4r3+Ldrxek9+P4YjqzvoX8KVSqbdbNEiGCcHLx727p3+3vFx2L7duV6xIr6xpbJQkumVwDuBi3HKPGzgWkRkSg2jL7GBh6igk3mblgNvDLuPY06VB7Wq8hAJWXAy3dcXeTKd3evMTGdXKZl2s4VBm3Hs2TP9mS3H7vobPxv6Fk9yHs9VXEpl5ZmJCTIFhZJMXw0stNaOxDsYEXGXha1P8Q4+AMDfdrydcJNpnw+OO+sPlUyLhMETtK17XxQbeuR7nWQ6r1Y1025WVeX/R5jX6/8709Hhf2wyffc+xLXcxbXcxT0FHwS+k9BYU0koyfQWoBydeigiYcqtcX7x5vWHXzTd1WkZHQUwlJZCUVHsYhNxuw95v0wRu/DQx8iur8CKhoj6KRh0kumC0zQz7WbG+Es9dmwZ4iz+Tt9Xnqdq7ijcdNMp9+ZtfGKi3X36KxMZZsoJJZmuBXYYYzbi1Exba+1V8QtLRNygMOgXb+FQ5zR3Tq5zUwuDLOc4teweXwXcG8PoRNzt0t7fsZxnANh48GP4d7gNX/Gw8w/hIiXTrrd4MQxsOcDjnA9fxf+R4Ec+cnK9h89HbctTE5c5Gb44PJRkOnhHHQOcD7wtPuGIiJsUz3VmpotHwp+Z7tl1nCUMM58DDJrqWIYm4nrD+R4IlHeMdERe5/F4bhMVo4uooJOaBZUxik5S1fLl8HsW0YuHUvr8tXb790NDg3PTjh0Uj/YA0Eo19Rdl5smHJ4RyAuIjxpjVwNuBtwD7gO/FOzARSX+l851kunQs/Jnpgb3O6sN+j04/FAnHaIFTND3WFXkyfSM/ZiDQ7l0QZVCS8pYvBx/ZPMW5XMb9/gefeOKkZHrs0ScmEsgnOY+LVk2xSjFDTHnuozFmqTHm88aY7cBtwEHAWGvXW2u/nbAIRSRtlTU4yXS57WJ8zIb1+pGDzurDkVlafSgSjrFCJ5kejzCZHhkxDAQy6exs/+I0cbdly/x/PklQ6caTT550T9/vH5hob684P+KdYtxiukPUdwCXAG+w1p4fSKDHExOWiLhBdlE+/fhXDeYwTs+h8H6h26POzLSvWjPTIuHwFTnJtK83smTa682daFdUTL1NmrhHo/8ARJ4gaFHh44877fFxiv7214nLtrWvSVBkqWu6ZPoa4BjwsDHmB8aYS/DXTIuIhKwn21mw1Ls/vLrprDZnZjqrTjPTIuGwJUF740WYTPf2OtWgM5wsLS5RXAzz58PTnMP4iTRx82Zn0/+NG8kf8P8sP0IdVetXJinS1DFlMm2t/b219jpgGdAMfAyoNcZ81xhzWYLiE5E015/n/Ab2HgwvmS7odmam8+drZlokLDHYaDp3Zwv/j8/ySb7KFVn3xSgwSXWnnw49lPMYF/gfsBb+8Ad/++67J+67j9ey9izNs043Mw2AtbbfWvsLa+3rgbnAZuAzcY9MRFxhoMCZmR44FN4ixGKvMzNdtEAz0yLhMGVOMp01EFkyXbJnL5/lS3yVT/HGnp/GKDJJdWvX+v+8h6BdkANJtH30Mech3siaNYmMLDXNmEwHs9Z2Wmu/b63VUeIiEpJdc9bza67jO3yQtuzwEuKyYSeZLm9UMi0SjuygZDo7wmTadHsn2qMl2mM6U5x1lv/Pk5Lp++6DvXvZcttjNPEwt/Ehtp52OTU1yYkxlYSyz3TMGWO+CrwBGAH2AO+x1nYnIxYRia+HLriF773gb98WxgmG4+NQPe6UeVSsUJmHSDhyZjnJdO5gZMl0Vo/zuvEyFU1nihPJ9D4W8lD2pVw8/gAUFsL+/Tzy/EIeoYlHaOK685MbZ6oIa2Y6hu4HzrDWngnsAj6bpDhEJM4qgiazOsOo8ujY76WEfgCGySOvuizGkYm4m2/FGfwrn+eTfJV7qt4bUR+5wQsXZ2lmOlPMncvEjPOt4//C8MJl8MwzsH49jz7q3HfhhcmJL9UkJZm21v7VWjsWuHwKfy22iLhQpMn00d5iajnGmWxhw/y/ak8ukTBlrVjGrfwrX+eT3J/3uoj6yO/vdfqrUjKdKYyBV7zC336MC/nvf3oeVqzAWk5Kpi+6KDnxpZpkzUwHuxH4c7KDEJH4iDSZPnbc0EotWzmTQ4v0E1skXMGbeXi9U983nfxBJ5nOqVEynUnWr3faDz7q329882Zob/c/VlnpPy1R4lgzbYx5AJisyPFma+09gXtuBsaAX0zTzwZgA0B1dTXNzc2xD1ZSgtfr1fi60NBzPXyJJ5hFF8OPzKW5ObTPBZuba4ETP6mP09y8PW4xSnT03k1NbW15EDh4o6NjmObmJ6d/wSTKB53lTIcGujTOLjPde7e0tARYB8Cf/jTCgw8+yR13zAcaAFiz5jiPPqqfyxDHZNpae+l0zxtj3g28HrjEWjvlGcPW2tuB2wEaGxttU1NTLMOUFNLc3IzG130qn/0LK/kKAM92X8K6ps+H9LpnnnHaq1bV0tSk3TxSld67qanXmVRmeDg/ojHaM9oz0V554TpWN50Rg8gkVUz33r3wQrj1VjhyBHp68hgdvYgnnnCef//79XP5hKSUeRhjXgt8GrjSWjuQjBhEJDGK5jofDRcNhV7n4d3bSi3HyGaM2drIQyRsxYU+7uFKHmI9D/WfjW98ynmrKXnGnJnpknkq88gkWVnwlrc415dfDnv2+Nvl5f5r8UtWzfRtgAe43xiz2RjzvSTFISJxVlLvbKdVMhr6CYgXPngLx6hjhDwu2PGDeIQm4mrZuVlcxl9ZTzNns5H+9sHwOrCWcp/znvXM09Z4meZ975t87feNN/qPHRe/pOwzba1dnIzvKyKJV9rgzGaV+TqxNrSNOfJ7/Ae2ZGEprCuPV3girtZvPBTYYX/7WB+e2tA3e/eNjvMffJxZdFFKL9fMKYxXmJKiVqyA97wHfvxj57GaGvisNjQ+SVKSaRHJHIWznf2hy+jF2ztOSVn2jK8r8ToHtpQsUl2eSCQGsj1Ujvm3Xxhs7QNCfy/1DebwWb4E+Gch35YXjwgl1d12G1gLd9wBixfDz34GVVXJjiq1pMLWeCLiZtnZdBtnZrm7JbTDTsuDjxJfpqJpkUgM5JROtP3JdOi6gqqyKlQunbEKC/0z06OjsHMnnHNOsiNKPUqmRSTuerOdZLpv/8yLEEdGoMbnzEzPWqaZaZFIDOU5m02PdISXTAfvC69kWrKUMU5J/2tEJO76cp1kuv/QzIsQW/d6Kca/0c8Q+WTPKp3hFSIymZH8oGS6vXeaO08VPDM9S2sPRaakmmkRibuB/FIIbCQwdGTmmenOF48x90Q7dzZzdJS4SERGC5xkeqwrvJnpvAf/zK/5KZ1U0Df4auBNMY5OxB2UTItI3A0WOIsQh4/NPDPdu8sp8egpqmNOXKIScb/xQieZHu8OL5nO3bGV67gTgL/0F6NkWmRySqZFJO62NlzME8dW0EY1S4pXc8kM9w/uc5LpwTItPhSJ1HiRk0z7esJLpm1Q0bRvloqmRaaiZFpE4m7b6jfwvacWAfDRECo2xg4enWiPViqZFomULXGSadsbXjKd1eUk06ZSybTIVJRMi0jclZWNTLTb22e+f7BriEEKKGQIO7sujpGJuNuBtW/khr8tpA8P5yw6g6YwXpvT6yTT2VVKpkWmomRaROKuvHx0ot3WNvP9P6v9J67hk3jo447r4hiYiMv1N67lDtYCUJ0b3mtz+531DXmzlUyLTEXJtIjEXbjJ9LFjAIY+SqleFLewRFzP41R50BdelQeFg87MdEGd9sYTmYqSaRGJu9njh/kVX6WaNoZ3lAO/nfb+o07JNHWq8hCJWHAy7fWG99riYSeZLq7XzLTIVJRMi0jclZaO8UZ+A8DRwdlYC1NtHW3tiZlpv9lafygSsZISpx3uzHTpmJNMe+YrmRaZipJpEYm7rBrnN3oV7Xj7LJ7SybPprg4fl4/cyzFm4y2eTXFxQ4KiFHGf6o4dPM+1eOij/dmFwEOhvXB0lBLrn8oeI5vyeTqFVGQqSqZFJP7y8+gzHjy2j1zGOLK3G8/qyWswW19s526uBqB7cBYw84mJIjK5omLDErYBYIZCX4E4cryLvEC7m3Iqy3QKqchUlEyLSEL05FbjGfF/zty9u535UyTTPTucgumugjrKExKdiDsV1TpF04Xjodd5dA8X8jm+TwWdFBUZPq9cWmRKSqZFJCG8hVUwsheA/pY2YMnk9+12Cqa9JSqYFolGcDJd7As9me4c9fADNgCw5DT4fMwjE3GPrGQHICKZYbC42mkfnPrklpEDTjI9PEvJtEg0SmqLJ9rFDGDHxkN6XdBJ4lRo7aHItJRMi0hCjJZVOe0jU282PX7YSaZ91UqmRaKRm59FH84C4MG20PbH63LOa2GWtpgWmZaSaRFJCF+lMzPta506mc5uc5LprNO0ybRItPqznFKPgeOhlXpoZlokdEqmRSQhTI2TTGd1TF3mkd/pLEDMn6+ZaZFoDUSQTM/7/bd4hldwH6/h4s674hWaiCtoAaKIJERunVPmkds99cx0Sd+RiXbpUiXTItEazPXAWKDdGloyXXjoJV7BswD80b4uXqGJuEJSkmljzBeAqwAf0ArcYK09Mv2rRCSdZV14Phv+6/u0U0Vu6RIumeQea6F66NDEdeWa+sQFKOJSw3keGPS3Q02ms3qcOo/sKtV5iEwnWTPTX7XW/guAMeYm/LvufCBJsYhIAlScu5QfsBSAuu7J72k77uM0nGS6pHFuIkITcbWRglLo8beH20NLpvP6nGQ6p0bJtMh0klIzba3tDbosBmwy4hCRxJkdVLFx/DiMT7JD15GX+mmmie0s42huPRQXn3qTiITl3rX/yoU8who2sW/BxSG9Jn/ASaYL5iiZFplO0mqmjTH/tVfjgAAAE2NJREFUF3gX/n8vr09WHCKSGHl5UFUF7e3g80FrK9S9bLOO/Z0e3shfAXjtJfDnJMQp4jbdDat5LNBuHwntNcVDTjJdNFfJtMh04pZMG2MeACZbPXSztfYea+3NwM3GmM8CHwZumaKfDeA/hqm6uprm5uY4RSzJ5vV6Nb4udWJsS0vX0d5eAlj+cO+zLG3sP+m+hx46jRMnI2ZnH6G5eVfig5Ww6b2b2rq7FwLzANi8eS/NzQdmfM3qUSeZfqlzP73NWtbkRnrvxkbckmlr7aUh3vpL4I9MkUxba28HbgdobGy0TU1NMYlPUk9zczMaX3c6MbbfGb6RxTxCHUf5+8CjXPCy8b7vPqd9zjlzaGqak9hAJSJ676a2J5+EX/3K366sXEhT08LpXzA+js/nLGxYf/XF1NRlxzFCSRa9d2MjKTXTxpglQZdXAjuSEYeIJFYdR1nEXooYZGDP0VOeP3jQaddrIw+RmCgr8/+ZzRgDHYMz3u/r6iErsJSpi3JmVSmRFplOsmqmv2SMacS/Nd5+tJOHSEYYrayDw/72yP5Tk+llf/9vPkYrB6lnYfGFQG1iAxRxodN3/pZebsCDl8cfeBv+D4Sn5t3fQWmg3WUqmJUb9xBF0lpSkmlr7TXJ+L4iklx2dh0872+bo6fWYF7W8gPO4VEADnrvR8m0SPQKyvLx4AUgb6Bnxvv79ndOJNO92eVxjEzEHXQCoogkTE69s31HTtvJM9M+H1QNBx3Yslp1HiKxkF9d5rSHptjkPUhrxTLezf1U0El5dWDRkohMScm0iCRMyVJnQWFB18kz00cO+ZgbdGBL0ZLTEhaXiJsV1DrJdNHIzDPTbSNlPIh/D4G1dV1xi0vELZKyAFFEMlPV2nlO27sfG3Rc04Gnj5KPfxPcnpwKKClJdHgirlQ0xynVKBqbOZnudHbFw+MZjUdIIq6iZFpEEqZsVcNEe77dR0e7k013btwz0W4vW5zIsERcreQ0Z2baMz5zmUdwMl1aOhaPkERcRcm0iCSMqaqkP8s/4+zBy4HNzm/toW27J9oDdYsSHpuIW3nmePBh/G282NHpE2TNTIuER8m0iCSOMbQVN0xcdvy9xXlqrzMzzWLNTIvESm5+Fr0T+3OA90jvtPdf9qsbaKOKnSzllV33xzs8kbSnZFpEEspb1eC0t7VMtEuOOTPTRSs1My0SS94sp9Sj/8j0ddP5Pa1U0cFSXqKoUDPTIjNRMi0iCbXlmi9wFs9SSTt/KX4TANZCdZ8zM111rmamRWKpP8dJpgeOTF83XdDv1HmYSk/cYhJxCyXTIpJQngtWs4mz6KSSl3b76zgP7Lcs9Dkz06VrlUyLxNJAnrOjx9Dx6Wemi4acZDqrqihuMYm4hfaZFpGEWrHCaW/b5v9z6+ZxHuQWVrKVs8p2s6q2JjnBibjU11b9nEcfz6KHMv5nQQkrprm3ZNRJpnNri+MfnEiaUzItIgm1cCEUFsLgILS2+r+2vJDDN/kYAB95F/ynSXKQIi4zVDt/4kiknr5pbvT5KB13DmrJqy2Ma1wibqAyDxFJqKwsOP10qKaVJh7mxY39PPec8/yZZyYvNhG3KnNKpumdZjMP29NLNj7/fXgomRXnwERcQMm0iCTczw+tp5VaHuZiDvxhM3/7m/PcunXJi0vErUqdnfHomaZkeuCQU+LRaSrJz/fFMSoRd1AyLSIJl9tw2kT7wF0bOX7cfxJieTmsXJmsqETcq7J0lCraWMRuxvcfmvK+vv1OMt2bU5GI0ETSnpJpEUm4sstfOdG+pOM3tFHN3VzFF0+7jezsJAYm4lKv3PUT2qhhN0s4//5bprzPe8BJpgfyVOMhEgol0yKScJVXXzjRPo+nqKKDq7iXq8Z+m8SoRNwrt9rZGi+7f+o6j6FDbRPtgaKquMYk4hZKpkUk4czpK+j31J7yeMX1VyQhGhH3y692ViDmDkydTG9ffg0L2MvZPM29q/4lEaGJpD0l0yKSeFlZ5G648aSHfCaL/Pden6SARNytsM6ZmS4YmvoExNbeAlpYwEbOZmjR6YkITSTtKZkWkaTI+/iHseXOL/isf/gg1NUlMSIR9yqudxYTeoY7pryvzanyoLo6nhGJuIcObRGR5JgzB7NxI/zoR7BoEdxwQ7IjEnEtT0PlRLt0vHPK+5RMi4RPybSIJM/ixfDFLyY7ChHXK28ox4chC0uZ7cGOjmFyT00Bxg4cppQSeimlulpHkYqEIqllHsaYTxpjrDFGS4ZFRETiJK8giy6cre6Ct8AL9rkH19NDOUMUMG9oV6LCE0lrSUumjTH1wKuBA8mKQUREJFP05DilHr0tk9dNlw776zzyGaF0vvaZFglFMmemvwF8CrBJjEFERCQj9OU6yfTAgUmS6dFRynz+nT58GCoW6wREkVAkJZk2xlwJHLbWbknG9xcREck0A4UVjJHNcWrwdgyf8ryvtX2i3UElVbU6jlQkFHFbgGiMeQCYPclTNwP/DFwWYj8bgA0A1dXVNDc3xypESTFer1fj61IaW3fT+KaHLyz7Hg89MRcw3Nq/jZ6XjZnd0sL6QLvDVPPCk80aW5fT+MZG3JJpa+2lkz1ujFkJLAC2GGMA5gKbjDFnW2uPTdLP7cDtAI2NjbapqSleIUuSNTc3o/F1J42tu2l808PPGoEn/O3TTjuDlw/ZwZYHJ9q9+dU0NTVpbF1O4xsbCd8az1q7Fag5cW2MaQHWWWvbp3yRiIiIRKUiqAS6c5LNPPpbnE2mvUXaZFokVDoBUUREJAMEH8ISfDjLCcOHnAeHPUqmRUKV9ENbrLUNyY5BRETE7epK+jiLnVTSQdnWYuD8k54fO+Yk06OzlEyLhCrpybSIiIjE35KOp3g2sPZ/y8b1wEMn3xA0XW0rlUyLhEplHiIiIhmgpME5bLhkoPWU58e8w4zh3w4vq642YXGJpDsl0yIiIhmgdKmzW235yPFTnv/Wqh+TzzB1HKHvgisSGZpIWlMyLSIikgEqGqvxYQCo9LVjR0ZPev7IEfCRzTHqqFlYkowQRdKSkmkREZEMUFKeQxtOLfRAy8mlHkeOOO05cxIVlUj6UzItIiKSAYyBjlyn1KN7x8nnpB096rTr6hIVlUj6UzItIiKSIXoKnGTau9tJpr0t7bzKex9nsoU5+R2UlycjOpH0pK3xREREMkS/pxb6/O2h/c4ixO6/PM19vB6AR7Muw5i/JCM8kbSkmWkREZEMMVzuzEyPHnJmpvt3OzUeXo9qPETCoWRaREQkQ/hqnGTaHnWS6eF9zurD4QqtPhQJh8o8REREMkRWwzz2sJDj1HLUzuMVgcd9h51k2lermWmRcCiZFhERyRDDb3gzi3/6ZgBeVwnXBB4vONoycY+ZPy/hcYmkM5V5iIiIZIj6eqd98KDT9nS2TLTzGhckLiARF1AyLSIikiEmTaatpdK7f+LxqrPmJzYokTSnZFpERCRD1NRAbq6/3dUF/f3AsWMU2CEAOplF/RllyQtQJA2pZlpERCRDZGXB66qepuboZhpo4fjD76S2oJviwPMtZgGrtf5QJCxKpkVERDLIh8e+ySX8GoBtjzTSVWonkuljxYvI0mfWImHRW0ZERCSDeGsWTrTHX9zF4f5yHmI9R5lNR82KJEYmkp6UTIuIiGSQkSWnT7Rzd23jwdKruYSHmMNRNr3uX5IYmUh6UjItIiKSQUrOPWOiXXFkGy+84Dy3/IzsJEQkkt6UTIuIiGSQuZc0MoY/aZ49sI/9W3snnluhKg+RsCmZFhERySBLzshnKysnriu2Nk+0lUyLhC8pybQx5lZjzGFjzObA1xXJiENERCTTFBTApsrLJq7v5Squ5+ecu7idiookBiaSppI5M/0Na+3qwNefkhiHiIhIRhm46PKTrn/Ou7jqjD1JikYkvanMQ0REJMMsePeFbGTdxPXDNLHwrWcnMSKR9JXMQ1s+bIx5F/As8AlrbVcSYxEREckYl16WxXmz/8ANx77IEAX8qPZmnr/KJDsskbRkrLXx6diYB4DZkzx1M/AU0A5Y4AtAnbX2xin62QBsAKiurj7rzjvvjEu8knxer5eSkpJkhyFxoLF1N41venrxxVK+/vWlZGdbPvGJXTQ29p1yj8bW3TS+01u/fv3frbXrZrovbsl0qIwxDcD/WmvPmOFWGhsb7c6dO+MekyRHc3MzTU1NyQ5D4kBj624aX/fS2Lqbxnd6xpiQkulk7eZRF3R5NbAtGXGIiIiIiEQjWTXTXzHGrMZf5tECvD9JcYiIiIiIRCwpybS19p3J+L4iIiIiIrGkrfFERERERCKkZFpEREREJEJKpkVEREREIqRkWkREREQkQkqmRUREREQipGRaRERERCRCSqZFRERERCKU9OPEw2GM6QN0nrh7VQHtyQ5C4kJj624aX/fS2Lqbxnd686211TPdlKwTECO1M5Qz0iU9GWOe1fi6k8bW3TS+7qWxdTeNb2yozENEREREJEJKpkVEREREIpRuyfTtyQ5A4krj614aW3fT+LqXxtbdNL4xkFYLEEVEREREUkm6zUyLiIiIiKSMlEymjTGvNcbsNMbsNsZ8ZpLn840xvwk8/7QxpiHxUUokQhjbjxtjXjTGPG+MedAYMz8ZcUpkZhrfoPvebIyxxhitIk8ToYytMeYtgffvC8aYXyY6RolcCD+b5xljHjbGPBf4+XxFMuKU8BljfmyMaTXGbJvieWOM+c/A2D9vjFmb6BjTXcol08aYbOC/gMuBFcDbjDErXnbbe4Eua+1i4BvAlxMbpUQixLF9DlhnrT0TuAv4SmKjlEiFOL4YYzzATcDTiY1QIhXK2BpjlgCfBV5lrT0d+MeEByoRCfG9+zngTmvtGuCtwHcSG6VE4afAa6d5/nJgSeBrA/DdBMTkKimXTANnA7uttXuttSPAr4GrXnbPVcAdgfZdwCXGGJPAGCUyM46ttfZha+1A4PIpYG6CY5TIhfLeBfgC/n8kDSUyOIlKKGP7PuC/rLVdANba1gTHKJELZXwtUBpolwFHEhifRMFa+yjQOc0tVwE/s35PAeXGmLrEROcOqZhMnwYcDLo+FHhs0nustWNAD1CZkOgkGqGMbbD3An+Oa0QSSzOOrzFmDVBvrf3fRAYmUQvlvbsUWGqMedwY85QxZrqZMEktoYzvrcD1xphDwJ+AjyQmNEmAcH83y8uk4gmIk80wv3zLkVDukdQT8rgZY64H1gEXxTUiiaVpx9cYk4W/LOuGRAUkMRPKezcH/8fETfg/UXrMGHOGtbY7zrFJ9EIZ37cBP7XWft0Ycx7w88D4+uIfnsSZcqoopeLM9CGgPuh6Lqd+nDRxjzEmB/9HTtN9hCGpIZSxxRhzKXAzcKW1djhBsUn0ZhpfD3AG0GyMaQHOBe7VIsS0EOrP5XustaPW2n3ATvzJtaS+UMb3vcCdANbaJ4ECoCoh0Um8hfS7WaaWisn0RmCJMWaBMSYP/0KHe192z73AuwPtNwMPWW2YnQ5mHNtAGcD38SfSqrlML9OOr7W2x1pbZa1tsNY24K+Jv9Ja+2xywpUwhPJz+W5gPYAxpgp/2cfehEYpkQplfA8AlwAYY5bjT6bbEhqlxMu9wLsCu3qcC/RYa48mO6h0knJlHtbaMWPMh4G/ANnAj621Lxhj/g141lp7L/Aj/B8x7cY/I/3W5EUsoQpxbL8KlAD/E1hTesBae2XSgpaQhTi+koZCHNu/AJcZY14ExoF/stZ2JC9qCVWI4/sJ4AfGmI/hLwG4QZNY6cEY8yv85VdVgZr3W4BcAGvt9/DXwF8B7AYGgPckJ9L0pRMQRUREREQilIplHiIiIiIiaUHJtIiIiIhIhJRMi4iIiIhESMm0iIiIiEiElEyLiIiIiERIybSIiIiISISUTIuIpBBjTKUxZnPg65gx5nDQ9RNx+p5rjDE/nOb5amPMffH43iIi6S7lDm0REclkgYNOVgMYY24FvNbar8X52/4z8O/TxNRmjDlqjHmVtfbxOMciIpJWNDMtIpImjDHewJ9NxphHjDF3GmN2GWO+ZIx5hzHmGWPMVmPMosB91caY3xpjNga+XjVJnx7gTGvtlsD1RUEz4c8Fngf/ceHvSNB/qohI2lAyLSKSnlYBHwVWAu8EllprzwZ+CHwkcM+3gG9Ya18BXBN47uXWAduCrj8JfMhauxq4ABgMPP5s4FpERIKozENEJD1ttNYeBTDG7AH+Gnh8K7A+0L4UWGGMOfGaUmOMx1rbF9RPHdAWdP048B/GmF8Av7PWHgo83grMif1/hohIelMyLSKSnoaD2r6gax/Oz/Ys4Dxr7SBTGwQKTlxYa79kjPkjcAXwlDHmUmvtjsA90/UjIpKRVOYhIuJefwU+fOLCGLN6knu2A4uD7llkrd1qrf0y/tKOZYGnlnJyOYiIiKBkWkTEzW4C1hljnjfGvAh84OU3BGady4IWGv6jMWabMWYL/pnoPwceXw/8MRFBi4ikE2OtTXYMIiKSRMaYjwF91trp9pp+FLjKWtuVuMhERFKfZqZFROS7nFyDfRJjTDXwH0qkRUROpZlpEREREZEIaWZaRERERCRCSqZFRERERCKkZFpEREREJEJKpkVEREREIqRkWkREREQkQv8fQ6ikaP0wFzYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Nlam = 16\n",
    "dx = vs0 / (Nlam * 2 * f0)\n",
    "dz   = dx    # grid point distance in z-direction (m)\n",
    "\n",
    "op = 2       # order of spatial FD operator\n",
    "\n",
    "# calculate dt according to CFL criterion\n",
    "dt = dx / (np.sqrt(2) * vs0) \n",
    "\n",
    "vy = FD_2D_SH_JIT(dt,dx,dz,f0,xsrc,zsrc,op)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... which leads to an improved waveform fit of the boundary reflections."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What we learned:\n",
    "\n",
    "* Estimating an analytical solution to the quarter-plane problem using the image method and superposition principle of linear partial differential equations\n",
    "* To accurately model the boundary reflections, we have to use at least $N_\\lambda = 16$ gridpoints per wavelength for the 2nd order FD operator."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
